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[57] ABSTRACT 

An embodiment of the present invention relates to a 
worldwide network of differential GPS reference sta- 
tions (NDGPS) that continually track the entire GPS 
satellite constellation and provide interpolations of ref- 
erence station corrections tailored for particular user 
locations between the reference stations Each reference 
station takes real-time ionospheric measurements with 
codeless cross-correlating dual-frequency carrier GPS 
receivers and computes real-time orbit ephemerides 
independently. An absolute pseudorange correction 
(PRC) is defined for each satellite as a function of a 
particular user’s location. A map of the function is con- 
structed, with “iso-PRC” contours. The network mea- 
sures the PRCs at a few points, so-called reference 
stations and constructs an iso-PRC map for each satel- 
lite. Corrections are interpolated for each user’s site on 
a subscription basis. The data bandwidths are kept to a 
minimum by transmitting information that cannot be 
obtained directly by the user and by updating informa- 
tion by classes and according to how quickly each class 
of data goes stale given the realities of the GPS system. 
Sub-decimeter-level kinematic accuracy over a given 
area is accomplished by establishing a mini-fiducial 
network. 

35 Claims, 6 Drawing Sheets 
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NETWORKED DIFFERENTIAL GPS SYSTEM 

This invention was made with Government support 
under contract NAS9-18360 awarded by NASA. The 5 
Government has certain rights in this invention. 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The invention relates generally to the global position- 10 
ing system (GPS) and specifically to methods and de- 
vices that allow a user to reduce GPS navigation errors. 

2. Description of the Prior Art 

The global positioning system (GPS) is a United 
States funded satellite system consisting of sixteen to 15 
twenty-four satellites in a constellation that beams 
highly accurate timed signals to earth. GPS receivers 
can process these signals from anywhere in the world 
and provide a user with position information that can be 
accurate to within twenty meters. Each GPS satellite 20 
moves through a known orbit and the time that individ- 
ual GPS satellite transmissions take to reach a user’s 
position are used in solving a triangulation problem. 
The individual signal paths from each satellite to the 
user’s position must pass through the troposphere and 2 
ionosphere. Position solution errors can occur when 
there is not accurately accounted for radio signal propa- 
gation delays caused by the troposphere and iono- 
sphere. The ionosphere introduces its maximum delays 
about 2:00 PM local time and has practically no effect at 
night. Thus, the local time may be important to a posi- 
tion solution. The troposphere can introduce delays that 
are a function of barometric pressure, temperature, 
humidity and other weather variables. Variations in the 35 
actual orbit of each GPS satellite can also have an im- 
pact on position solution accuracy. Drifts in the clock of 
each GPS satellite will also impact user position solu- 
tion accuracy. Service bureaus and agencies have been 
established to sell or otherwise provide differential data. ^ 
Such data range from real-time local information to 
sometimes very complex mathematical models based on 
long-term observations. US Government agencies issue 
orbit correction information on a satellite-by-satellite 
basis. 45 

Dual frequency carrier GPS receivers continuously 
track P-eode LI and L2 carriers of a GPS satellite to 
generate accumulated delta-range measurements 
(ADR) and at the same time track LI C/A-code to 
generate code phase measurements. Each carrier is 50 
modulated with codes that leave the GPS satellite at the 
same clock time. Since the ionosphere produces differ- 
ent delays for radio carriers passing through it having 
different radio frequencies, dual carrier receivers can be 
used to obtain real-time measurements of ionospheric 55 
delays at a user’s particular position. (LI is typically 
1575.42 MHz and L2 is typically 1227.6 MHz.) The LI 
and L2 ADR measurements are combined to generate a 
new LI ADR measurement that has an ionospheric 
delay of the same sign as the ionospheric delay in the LI 60 
pseudorange. Accurate ionospheric delay figures, if 
used in a position solution, can help produce much 
better position solutions. Without such real-time iono- 
spheric delay measurements, mathematical models or 
measurements taken by third parties (which can be old) 65 
must be used instead. The communication of this infor- 
mation to a user’s site can be costly and require wide 
communications channel bandwidths. 
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Commercial and private users have been able to make 
use of GPS satellite transmissions even though they use 
so-called “unauthorized” receivers. “Authorized” re- 
ceivers are able to receive special information (P-code) 
that can remarkably improve position solution accu- 
racy. Since an enemy could use GPS to fly and target 
missiles and artillery shells to better than one meter of 
accuracy, the United States Department of Defense 
occasionally engages a selective availability (SA) mode, 
spoofing and encryption of the GPS signals that deliber- 
ately introduce solution errors into user receivers that 
lack an authorization code needed to see through the 
masking. A “codeless” unauthorized GPS receiver is 
able to cross-correlate the LI and L2 signals and extract 
enough information to measure the ionospheric delay 
(although not enough to overcome SA or to read P- 
code encrypted by Y-code to secure transmissions.) 

With differential GPS, a stationary reference receiver 
is placed at a very accurately known point location. The 
reference receiver generates corrections which are sent 
to a transmitter, which in turn broadcasts the correc- 
tions to users within the area of the transmission broad- 
cast. A differential GPS user receives these corrections 
through a radio/modem and applies them to the direct 
GPS measurements. This gives the user a position mea- 
surement of very high accuracy, e.g., from one meter to 
ten meters. 

However, as a user’s receiver is moved away from a 
reference receiver, the corrected accuracy will gradu- 
ally deteriorate. Thus, the region of usefulness is finite. 
By employing a network of reference stations, the accu- 
racy achievable by users near a reference station can be 
had throughout a region of several hundred miles. 

While advantages of differential GPS networks have 
been well known, the basic problems to realization of 
them are many. A means of reducing the communica- 
tions volume between the reference stations is needed in 
order to keep the cost of service down and thus to be 
competitive. A means of reducing the communications 
volume between the reference station and the user pop- 
ulation is needed. Also needed are effective techniques 
and algorithms that enable a minimum number of refer- 
ence stations to be employed. Effective techniques and 
algorithms are needed that provide high accuracy 
throughout a large operating region, for example, 
worldwide in a band of latitudes that include most of 
the populated areas. 

A paper written by Alison Brown, titled “Extended 
Differential GPS,” published by Navigation: Journal of 
the Institute of Navigation, in Vol. 36, No. 3, Fall 1989, 
pp. 265-285, presents an extended differential GPS 
concept that uses a network of differential GPS stations. 
A differential GPS message can be computed that pro- 
vides corrections for the different components of the 
pseudorange error. This is reported to extend differen- 
tial GPS navigation accuracy to one thousand nautical 
miles from a master differential station. 

NETWORKED DIFFERENTIAL GPS 

GPS navigation error results from a variety of 
sources. In general, the worst error contributors are 
common, virtually constant and exist over large areas. 
Included are satellite orbit and satellite clock errors 
(especially that due to the intentional degradation of 
selective availability), ionospheric error and to a lesser 
extent, tropospheric error. Differential GPS (DG PS) is 
a prior art method of eliminating such errors. It works 
by measuring the GPS signals with a receiver at a 
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known location (a reference station). The actual posi- 
tion is compared to the position solution. A difference in 
the positions can be due to ionospheric delays, satellite 
clock bias and orbit ephemerides. The net difference of 
all these is broadcast via any of several means to a local 5 
group of GPS receivers for correction. Differential 
GPS works quite well to produce meter-level accurate 
navigation as far as 100 kilometers from the reference 
station. The 100 kilometers limit is roughly consistent 
with line-of-sight high frequency (HF) differential data 10 
links, which are convenient for the high data rates (e.g., 
300 bits per second) desirable for differential GPS. Con- 
sequently, simple DGPS has always been satisfactory 
for a wide range of navigation and positioning applica- 
tions. However, the accuracy degrades as the distance 15 
from the reference station increases in amounts that are 
unacceptable to various users. In addition, many of the 
larger GPS errors will vary significantly over just a few 
hundred kilometers. It has been estimated that the 
DGPS degradation rate is approximately one centime- 20 
ter of accuracy per kilometer (km.) of distance from the 
reference station. 

There are at least three situations in which a simple 
DGPS solution is not sufficient. First, if the purpose is 
to provide DGPS service over a large area, establishing 25 
differential reference stations every 200 kilometers is 
extremely inefficient in terms of hardware and data link 
frequency allocation. Second, it is sometimes impossible 
to place a differential reference station within 100 kilo- 
meters of a working location. For example, oil explora- 30 
tion operations in the Gulf of Mexico and the North Sea 
need high levels of DGPS accuracy far from shore 
where the reference stations can be positioned. Also, 
there are certain users, such as the hydrographic sur- 
veyors servicing shipping channel dredges, that require 35 
decimeter level accuracy over a few dozen kilometers. 
Such levels of accuracy cannot be met by DGPS alone. 

An embodiment of the present invention, a net- 
worked DGPS (NDGPS) answers the above needs. 
NDGPS reference stations can be thousands of kilome- 40 
ters apart and still supply much better coverage and 
accuracy than DGPS stations spaced 200 kilometers 
apart. Long distance off-shore service is handled by 
encircling a work area with DGPS reference stations. 
Corrections are interpolated for each user’s site on a 45 
subscription basis. The data bandwidths are kept to a 
minimum by transmitting information that cannot be 
obtained directly by the user and by updating informa- 
tion by classes and according to how quickly each class 
of data goes stale given the realities of the GPS system. 50 
To provide sub-decimeter-level kinematic accuracy 
over a given area, a mini-fiducial (fixed reference) net- 
work can be established. 

GPS ERROR SOURCES AND SPATIAL 55 
DECORRELATION 

In the present invention, an absolute pseudorange 
correction (PRC) is defined for each satellite as a func- 
tion of a particular user’s location. A map of this func- 
tion is constructed, with “iso-PRC” contours much like 60 
isobar contours on a weather map. The object of the 
network is to measure the PRCs at a few points, so- 
called reference nodes and to construct an iso-PRC map 
for each satellite. Due to the temporal changes in the 
ionosphere and the continuous change in the orientation 65 
of the line-of-sight vector to the space vehicle (SV), this 
map must be frequently update.. The spacing between 
the reference nodes is kept as wide as possible, so that 
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any inferences about the character of the PRC function 
can be taken into account. If, for example, the PRC 
function varied unpredictably over a hundred kilome- 
ters, the reference nodes would have to be closely 
spaced. Fortunately, in general, this is not the case. 

The error sources unique to each GPS receiver are 
multipath, code noise and inter-channel bias (where 
applicable). These errors are uncorrelated between re- 
ceivers. The other error sources are correlated to some 
degree. DGPS provides pseudorange corrections de- 
signed to exactly compensate these errors at the GPS 
reference stations. However, as the user gets farther 
away from the reference station, the corrections be- 
come less exact. For each error there is some mecha- 
nism that controls the spatial decorrelation effect with 
distance from the reference station. Table I summarizes 
some of the major error sources and their spatial and 
temporal variations. The temporal decorrelation deter- 
mines the frequency at which the DGPS pseudorange 
corrections have to be updated. The spatially correlated 
GPS error are described herein. The spatial decorrela- 
tions determine the accuracy of DGPS correction that 
can be achieved. 

Satellite Clock Error — GPS navigation is basically 
hyperbolic navigation and relies on measuring pseudo- 
ranges from four or more transmitters (on satellites) that 
have known (albeit dynamic) positions and synchro- 
nized signal transmission. Small errors in the signal 
synchronization do occur and are calibrated for and 
broadcast. Under selective availability mode, the syn- 
chronization between the satellites is purposely made to 
wander. This effect, called clock dither and all other 
satellite synchronization errors, are completely corre- 
lated among all receivers that are tracking that satellite. 
Thus, the DGPS correction error is independent of the 
distance between the reference station and the user. 


TABLE I 
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Location 

(No Correlation) 

(No Correlation) 


Satellite Orbit Error — The orbits of the satellites are 
measured by U.S. Government agencies and broadcast 
subject to selective availability. The orbit parameters 
are purposefully misreported a small amount to cause a 
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controlled navigation error for the user, the so-called 
epsilon error. These slightly incorrect orbit parameters 
will cause both the reference station and user to use an 
incorrect satellite transmitter position in their computa- 
tions of a user position solution. Although the user and 5 
reference will have identical errors in computed satel- 
lite position, they will have slightly different errors in 
their respective computed ranges because of the differ- 
ences in their respective viewing angles to the GPS 
satellite. As the separation distances between reference 10 
and user grow, so will the differences in viewing angles 
and the differences between the computed range errors. 
The navigation error caused by satellite position error is 
highly correlated between viewers and very linear in 
nature (e.g., the error is roughly proportional to the 15 
distance from the reference station). 

Ionospheric Error — As the GPS signal travels 
through the ionosphere, it experiences a delay for 
which there must be compensation. This delay can 
either be measured with a receiver that is capable of 20 
dual frequency code measurement, or it can be mod- 
eled. For ionospheric model errors, the correlation 
distances mainly depend on the particular model being 
used. The ionospheric models that serve stand-alone ^ 
navigation, such as the GPS broadcast models of 1 
Klobuchar, Bent and the most recent version of the 
International Reference Ionosphere, generally use 
monthly averages for prediction and are believed to 
help remove 50% to 80% of the ionospheric delay un- 3Q 
certainty. However, relying solely on ionospheric mod- 
els can leave an error of one-third to one-half of the 
ionospheric delay. 

The following general observations can be made 
about ionospheric models: all models predict quiet iono- 35 
spheric conditions only, the mid-latitude trough which 
exhibits large horizontal gradients in electron density is 
not incorporated in these models, the models work well 
only for latitudes between ±20 degrees and ±60 de- 
grees, and they are poor predictors for equatorial and 4 0 
high-latitude regions. Since these very regions can be 
the principal areas of operation for certain users, these 
users are disaffected by the models. 

Tropospheric Error— The tropospheric delay in a 
GPS measurement is an unwanted delay in the code and 45 
carrier data which is introduced by the earth’s tropo- 
sphere, which extends from sea level to an altitude of 
approximately fifty kilometers. Since tropospheric 
delay is RF frequency independent, it can be modeled. 
The total delay consists of a dry and a wet component. 50 
The dry component, which accounts for about 90% of 
the total delay, can be reasonably well modeled without 
any meteorological data. The wet component, on the 
other hand, requires measurements of the local weather 
conditions along the lines-of-sight to the GPS satellites 55 
being tracked for maximum accuracy. The best models 
extrapolate surface weather measurements, such as 
pressure, temperature and humidity for the higher alti- 
tudes using standard atmospheric profiles. All tropo- 
spheric models are good above fifteen degree mask 60 
angles. Comparisons of tropospheric delay models rela- 
tive to ray-trace tropospheric delay estimates show that 
global models can be used to predict tropospheric delay 
to accuracies of 10 cm to 20 cm for elevation angles 
between 10° and 20°. Direct surface refractometry is a 65 
second method for dealing with tropospheric delays. 
The delays are measured directly along the line-of- 
sight. 
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For navigation networks, with hundreds of kilome- 
ters between reference stations, tropospheric error is 
more or less uncorrelated and hence considered to be 
reference station unique. For small networks, especially 
mini-fiducials, tropospheric model error between the 
reference stations and the users may be correlated and 
the network solution can be used to help improve user 
accuracy by including this term in the correction mix. 

Datum and Clock Error — By selecting the position 
coordinates of the network reference stations, a net- 
work can effectively establish its own datums. These 
datums are related to existing datums, such as WGS84. 
Survey errors in the reference station coordinates will 
cause a correlated error among network users. Gener- 
ally, a DGPS user will be navigating in both position 
and time, using the reference station position and clock 
as an absolute reference. This is also true for NDGPS, 
except that the network clock and datum are used as the 
absolute reference. 

Another source of error is the latency of corrections 
due to constraints in transmitting out the information in 
real-time. This is a very important consideration in 
choosing data links and designing information flow. It 
may or may not be correlated among the network par- 
ticipants, depending on the particular DGPS algorithm 
implementation in the user receiver. 

SPATIAL DECORRELATION OF ORBIT ERROR 

Nominal Orbit Errors— Using just a few dual-fre- 
quency reference stations and an extensive satellite orbit 
and clock dynamics model, U.S. DOD agencies has 
managed to predict orbit and clock so well that user 
range accuracy of one to three meters is common. Orbit 
position errors are decomposed into three directions 
defined by the orbit itself: radial, along-track and cross- 
track. Radial error is error in the orbital radius. Along- 
track error is error along the satellite velocity vector. 
Cross-track error is error orthogonal to the radial and 
along-track errors. Radial error will translate directly 
into stand-alone navigation error without a significant 
affect on differential operations. For every meter of 
error in the radial direction, there can be nearly a meter 
of range error for every GPS receiver. For users di- 
rectly under a particular GPS satellite, there would be 
exactly a meter of range error. For GPS receivers that 
see the satellite on the horizon, about 0.97 meters of the 
error is observed due to viewing the radial error about 
15° off-axis. Consequently, the radial orbit error of a 
satellite position is viewed almost identically by all 
users. It is about equal to typical satellite clock bias 
error. This will not be true for along-track and cross- 
track errors. A meter of along-track or cross-track error 
(which are equivalent for terrestrial navigation pur- 
poses) is viewed by users 75° to 90° off-axis. Thus, it will 
not have much effect on stand-alone position accuracy. 

selective availability — U.S. DOD has reported that it 
will constrain the statistics of the horizontal error of a 
GPS-derived position to 100 meters 2 drms (roughly the 
95th percentile). This means that for 5% of the time the 
size of the position error is unconstrained. Further con- 
straints at higher quantile levels have been proposed, 
but there is currently no reasonable maximum limit in 
the size of GPS position error. 

There are two independent error sources that con- 
tribute to SA error, clock dither and orbit epsilon. 
Clock dither allows unpredictable clock synchroniza- 
tion wander of a few tens of nanoseconds. Thus, a sta- 
tionary receiver will find its GPS position wandering a 
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few tens of meters over a couple of minutes. Almost all 
of the selective availability experienced to date seems to 
have been clock dither. Clock dither is uniformly seen 
by all the GPS receivers, so DGPS can correct for it as 
well as can NDGPS. 5 

There are practical limits to epsilon and its spatial 
decorrelation. Assuming that clock dither and epsilon 
each contribute, there is about seventy meters of error 
at the 95th percentile level. In a simplest of assumptions, 
the orbit epsilon error is considered to be (statistically) 10 
equal in all directions. This can be validated by consid- 
ering the geometry of the satellite position error and 
how it translates into a pseudorange error. Many meters 
of along-track and cross-track orbit error are applied 
before achieving the assumed limit of seventy meters 15 
for the epsilon contribution to a horizontal error in 
stand-alone navigation. However, different users will 
view this error from significantly different angles, thus 
leading to a good level of spatial decorrelation. On the 
other hand, a radial error will introduce very little spa- 20 
tial decorrelation. Thus, compared to clock dither, the 
radial error is actually an inefficient method of selective 
availability implementation. 

A simple Monte Carlo simulation can illustrate the 
spatial decorrelation effects that are possible with epsi- 25 
Ion. If Gaussian errors with a standard deviation of 190 
meters are applied in both the cross-track and along 
track directions, with no error applied in the radial 
direction, a stand-alone horizontal navigation error of 
71 meters two drms results, which is roughly the level 30 
desired. Spatial decorrelations are measured by comput- 
ing the position error of a differentially corrected user 
one kilometer away from a reference station. For pur- 
poses of this example, the horizontal navigation error 
(after differential correction) is 1.3 cm two drms over a 35 
one kilometers baseline, or 13 ppm. The 95th percentile 
vertical error is 3.0 cm or thirty parts-per-million (ppm). 

At this level, epsilon becomes a significant error con- 
tributor in DGPS. For a separation distance of one 
hundred kilometers between a reference station and 40 
user, errors of 1.3 meters horizontal and three meters 
vertical 5% of the time can be expected if the broadcast 
orbit parameters with the above described epsilon are 
used in conjunction with NDGPS. Many survey opera- 
tions use post-processing, which allows orbit relaxation 45 
techniques or a replacement of the epsilon-con- 
taminated broadcast orbit parameters with precise orbit 
parameters obtained from some other source. Unfortu- 
nately, typical survey applications have a real-time need 
for accurate orbit parameters. High error levels are 50 
particularly disastrous for “instantaneous” carrier ambi- 
guity resolution methods, which are subject to false-fix 
and no-fix situations at high error levels. This makes 
ephemerides broadcast by various providers virtually 
useless in such applications because the ephemerides are 55 
stale. Code-averaging methods of carrier ambiguity 
resolution usually do not depend on the ephemeris to 
resolve the carrier ambiguity, but are still subject to 
spatial decorrelation error. 

Two choices exist for accurate orbit parameters in 60 
real-time surveys establish a mini-fiducial network, or 
estimate a set of orbit parameters, both can be practical 
solutions. The mini-fiducial network adds two more 
reference stations to an existing real-time survey setup 
and has the attendant equipment and data-link costs. 65 
Orbit parameter prediction requires a receiver tracking 
network and near real-time post-processing using exten- 
sive orbit models, so it is subject to error only in the 
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unlikely event that the satellite moves during the pre- 
diction interval. 

SPATIAL DECORRELATION OF IONOSPHERIC 
MODEL ERROR 

A single-frequency receiver relies on an ionospheric 
model to remove ionospheric delay from each pseudo- 
range. The actual delay depends on the density of iono- 
sphere encountered by the signal in transit and can 
throw off position solution accuracy many meters, espe- 
cially when the satellites being tracked are low on the 
horizon and cut through a maximum slice of the iono- 
sphere. A dual-frequency receiver can effectively mea- 
sure the delay making use of the second frequency. 
Dual-frequency receivers can be roughly categorized 
into three types: those that rely on knowledge of the 
P-code to measure the ionospheric delay; and those that 
do not measure the ionospheric delay at all, but rather 
track both carrier frequencies. Codeless (squaring) L2 
channel receivers may have long time constants and can 
have a lag error under active ionospheric conditions. 
Inter channel (L1/L2 channel) bias calibration may be 
necessary and must be possible. Conversion of dual-fre- 
quency measurements into ionospheric delay is not 
perfect, especially for satellites near the horizon. In 
general, the signal strength of P-code receivers is much 
higher and they track the ionosphere much more accu- 
rately and reliably, especially in active ionospheric con- 
ditions. 

With differential GPS, a measurement of the iono- 
sphere is made in the vicinity of the user. Depending on 
whether the receiver at the reference node is a single or 
dual-frequency type, the measurement is implicit in the 
pseudorange or can be separated out. An absolute iono- 
sphere model is not required. A model of ionospheric 
variation over the differential GPS coverage area can 
be enough. 

Ionosphere error is subject to occasionally very non- 
linear behavior, unlike range errors due to orbit parame- 
ter error, which are very linear over large distances. 
Various published results by expert researchers demon- 
strate that the range of ionospheric behavior in some 
measurements, such as those taken in Greenland, vary 
in decorrelation distance from 200 kilometers (under 
“disturbed” conditions to 1,000 kilometers (in the nor- 
mal case). Long-distance spatial decorrelation studies 
by Klobuchar indicate correlation distances in thou- 
sands of kilometers, with apparent differences observed 
in the east-west (3,000 kilometers) and north-south 
(1,800 kilometers) axes. Over moderate distances of a 
few hundred kilometers, the differential ionospheric 
variation has been measured in one instance to be on the 
order of 0.53 mm of differential vertical ionospheric 
delay per kilometer of separation distance, under day- 
time conditions. Nighttime conditions are usually more 
subdued. Other daytime results of 2.8 mm differential 
ionospheric delay error per kilometer of separation 
distance have also been reported in the literature. 

NETWORK ALGORITHMS 

From an algorithm standpoint, there are two basic 
varieties of networks. The simple type is the common- 
view network. The more complex type is the wide-area 
network. 

Common-View Network Algorithm Design— The 
common view network assumes that a differential GPS 
user will navigate using only satellites that are visible to 
all the reference stations. The network will not be very 
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large if the user is to be allowed track satellites very 
near the horizon. A mini-fiducial network covering an 
area 200 kilometers wide is a typical application of a 
common-view network. The elevation angle of a satel- 
lite near the horizon varies by roughly two degrees over 
the network. Since a common-view network area is 
small, it can be assumed that the PRCs are basically a 
linear function over the area. Across an area of such 
size, the linearity assumptions are reasonable for a sin- 
gle-frequency ionospheric model. It may extend to tro- 
pospheric models as well. The assumption may be over 
extended in the case of very low elevation satellites. 

The PRC applied by the common-view NDGPS user 
is a weighted average, or blend, of PRCs from the refer- 
ence stations. The weighting coefficients for a weight- 
ed-average are determined solely by the relative geome- 
try of the user and reference stations. The same 
weighting coefficients are used for all code and carrier 
data. If the user is located close to one of the reference 
stations the data from that node is more highly 
weighted. The weighting coefficients can either be de- 
termined analytically or statistically. In the latter case, a 
Gauss-Markov spatial decorrelation assumption is in- 
voked. Three weighting coefficients: ai, a2 and a3 are 
used for a three node network and are determined ana- 
lytically by solving the three equations: 

user latitude =2(a/) (node latitude)/; 

user longitude =2(a/) (node longitude); and 


l=Za i 

It is good practice to include more than three refer- 
ence stations in a network. Using a least-squares analysis 
techniques, the extra data is used for checking assump- 
tions about linearity and to check for reference station 
failures. The extra stations can also provide protection 
against equipment failure. 

In contrast to DGPS, which sends the PRC in a use- 
ful form to the user, the NDGPS user must compute its 
own blended PRC, because the weighting coefficients 
depend upon the user’s position. There are a couple of 
different methods to provide the user with the informa- 
tion to do this. One method, referred to as a decentral- 
ized broadcast method, has each of the three reference 
stations transmit normal, high rate DGPS streams. The 
user then computes the weighting coefficients aj, a2and 
a3 according to a rough estimate of the position and 
blends the three reference station PRCs before applying 
them to the pseudorange measurements. (The user can 
blend PRCs from more than three stations using a least- 
squares technique.) This method may be retrofitted to 
existing DGPS installations, because all the network- 
specific algorithms can be computed in the user’s GPS 
receiver. It is, however, very inefficient from the stand- 
point of frequency allocation, because three high-rate 
DGPS data streams are needed, which requires three 
times the communication bandwidth and so multiplies 
recurring operational costs. 

A more efficient method used in embodiments of the 
present invention, is a centralized broadcast method 
that integrates a network by transmitting a single 
NDGPS data stream. The single stream is a normal 
DGPS data stream with the PRCs for some arbitrary 
location and the north and east gradients of the PRC 
functions, relative to the arbitrary location. The DGPS 
data stream carries high-frequency information, which 
comprises mostly clock dither data. Typically, the in- 
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formation is carried by each of the DGPS data streams, 
resulting in a considerable waste of bandwidth. The 
gradient terms, due to orbit error and ionospheric error, 
vary relatively slowly and need only be updated every 
5 few minutes. The extra information added requires only 
a very low bandwidth. 

Two points need to be made about common- view 
systems that distinguish them from wide-area systems. 
First, a common-view system needs no synchronization 
10 between the reference node clocks. Whatever the indi- 
vidual reference clock errors are, they are averaged into 
pseudorange corrections to form a blended “network 
clock”, using the same weights that are used to blend 
the PRCs. The NDGPS user clock error is relative to 
15 the blended clock, just as the DGPS user clock error is 
relative to the DGPS reference station clock. Second, 
because the area of coverage is relatively small, there is 
no need to try to separate the satellite clock error from 
the satellite radial error (or more exactly, line-of-sight 
20 error). Consequently, the satellite position is not fully 
determined by the common-view network. 

In March of 1991, a series of tests were conducted by 
the inventors to evaluate the performance of a network 
25 of GPS reference stations as a source for differential 
D corrections. These corrections were used in a naviga- 
tion solution of a railroad utility vehicle. Four GPS 
reference stations surrounding a segment of the track 
were used in a Burlington Northern Railroad test on the 
Mesabi Iron Range in the northern Minnesota. The 
results obtained demonstrate that a DGPS network 
applied to a railroad environment can achieve the most 
stringent accuracy requirement, a resolution of parallel 
tracks a mere four meters apart. 

35 WIDE AREA NETWORK DESIGN 

Wide-area networks can be designed to be efficient in 
terms of the installed hardware and communications. 
They also provide a potential for complete observabil- 
4 q ity of the separate spatially correlated DGPS error 
sources. Wide-area networks typically span a much 
greater area than are covered by a single data link bea- 
con. Thus, satellite communication is a requirement for 
a wide-area system. A wide-area network is essentially 
45 a satellite tracking network. GPS adds unique capabili- 
ties to the tracking network. Each tracking station, or 
reference node, is relatively inexpensive, in terms of 
receiver hardware and each can track many satellites 
simultaneously. An abundance of data allows tracking 
50 stations to do time transfer, or reference node clock 
synchronization, at the same time they are tracking the 
satellites. Algorithmically, the simplest wide-area net- 
work is a patchwork of common-view networks. Each 
small network patch supports a NDGPS beacon that 
55 services a small area. No reference station clock syn- 
chronization, or any substantial integration of any kind, 
is necessary. At a station-to-station separation of 200 
kilometers, a network of less than 200 stations could 
cover the continental United States, providing universal 
60 common-view NDGPS PRCs of excellent accuracy for 
single-frequency users as well as PRC gradients for 
single-frequency survey applications. 

The next simplest approach algorithmically is a much 
sparser network on a much grander scale, the dual-fre- 
65 quency world-wide (WWDGPS) tracking network. A 
sample network with thirty-three reference stations 
located around the world generates enough measure- 
ments to instantaneously solve for satellite position, 
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satellite clock error and synchronize the reference sta- 
tions. In addition, there are enough redundant measure- 
ments to do validity testing and accommodate two or 
three reference station failures. WWDGPS navigation 
accuracy is estimated by the inventors to be on the 5 
order of one meter for a P-code system user. The same 
network gathers a sufficient number of absolute iono- 
spheric delay measurements to determine near real-time 
updates to a global ionospheric model (used by single- 
frequency users). The algorithm is a straightforward io 
application of least-squares, with 96 unknown satellite 
coordinates (24 satellites, each with three position er- 
rors and one clock error). Thirty-two reference clock 
errors are measured by over two hundred pseudo- 
ranges. An observation matrix for a least-squares esti- 15 
mator is created from the pseudorange equation, given 
range equation (SV I, reference node J): 
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errors (e.g., orbit and ionospheric errors) allows the 
NDGPS correction accuracy to be spatially indepen- 
dent relative to the reference point location. Since these 
specific errors vary slowly, lower frequency broadcasts 
of these corrections can be made to lower communica- 
tions costs. Wide-area networks also permit the use of a 
sparse matrix of reference stations. One of the tradi- 
tional disadvantages of wide-area networks are the 
higher communications costs since satellite data trans- 
missions may be required. 

Therefore, what is required for non-authorized users 
that require high accuracy GPS position determinations 
is a networked differential GPS system of reference and 
control stations that is capable of providing correction 
information regardless of a user’s particular location 
apart from a reference station and in spite of any selec- 
tive availability concerns. 


Rij= \SVPosi-RefPosj\ +SV Clocks Ref 
Clockj+ Noise 

The unknowns are three SV position errors for each 
SV, one SV clock error for each SV and one reference 
clock error for each reference station. (The equation 
assumes that the ionospheric and tropospheric delay 
errors are removed at each of the reference stations.) 
Satellite velocities and satellite reference frequencies 
are solved for with the same observation matrix using 
the carrier measurements. In a worldwide network with 
a little over thirty reference stations, all satellite position 
and clock errors can be determined directly from the 
instantaneous data. With fewer than twenty stations, 
there is not enough measurement data to resolve the 
system errors, so an orbital dynamics model is required 
to collect enough information for the estimator. In spite 
of the extra hardware cost of a larger, fully-measured 
system, there are several worthwhile advantages, in- 
cluding survivability in the face of reference station 
failures and complete visibility into most or all satellite 
soft failures, especially in the control systems. Even so, 
data-link loading from the reference stations to the cen- 
tral processing facility, which is a major recurring cost, 
is roughly the same. 

Once the satellites’ positions and clock errors have 
been estimated, satellite position and clock motion can 
be tracked continuously by a fully-measured network 
with dual-frequency (ionospherically corrected) carrier 
data that is commonly available in survey-quality re- 
ceivers without P-code capability. Consequently, a 
world-wide system can run accurately on dual-fre- 
quency carriers once satellite positions and clock esti- 
mates have been initialized. The difficulty is in finding 
these initial estimates, because they can not be easily 
calculated where direct ionospheric measurements are 
not available. The dual-carrier WWDGPS estimates 
converge fairly quickly using data collected for night- 
time ionospheric conditions. 

The same can hold true even for wide-area networks 
covering smaller portions of the globe. The difficulty 
with smaller networks is creating the initial estimate of 
the satellite position and satellite clock. This hurdle 
must be overcome more than once, because a regional 
system will not be able to track each satellite through its 
whole orbit and so each satellite returning after an ab- 
sence away from the network must be re-initialized. 

Networked DGPS provides local, regional, or 
worldwide differential corrections for use in real-time 
navigation or survey applications. With a wide-area 
network, identification of the key correlated DGPS 


SUMMARY OF THE PRESENT INVENTION 

It is therefore an object of the present invention to 
provide a system on a worldwide or regional basis that 
is capable of providing differential GPS correction 
information regardless of a user’s particular separation 
25 from any nearest reference station and in spite of any 
selective availability concerns. 

Briefly, an embodiment of the present invention is a 
network of differential GPS reference stations 
(NDGPS) that continually track the entire GPS satel- 
30 lite constellation and provide interpolations of reference 
station corrections tailored for particular user locations 
between the reference stations. Each reference station 
takes real-time ionospheric measurements with codeless 
cross-correlating dual-frequency carrier GPS receivers 
35 and computes real-time orbit ephemerides indepen- 
dently. An absolute pseudorange correction (PRC) is 
defined for each satellite as a function of a particular 
user’s location. A map of the function is constructed, 
with “iso-PRC” contours much like isobar contours on 
40 a weather map. The object of the network is to measure 
the PRCs at a few points, so-called reference stations 
and to construct an iso-PRC map for each satellite. 
Corrections are interpolated for each user’s site on a 
subscription basis. The data bandwidths are kept to a 
45 minimum by transmitting information that cannot be 
obtained directly by the user and by updating informa- 
tion by classes and according to how quickly each class 
of data goes stale given the realities of the GPS system. 
Sub-decimeter-level kinematic accuracy over a given 
area is accomplished by establishing a mini-fiducial 
network. 

Advantage of the present invention is that an 
NDGPS system is provided that is able to compensate 
55 for the problem of separation between a corresponding 
reference station and user. 

Another advantage of the present invention is that an 
NDGPS system is provided that eliminates the problem 
of DGPS users on the fringe seeing satellites not also 
60 visible to a single reference station. 

Another advantage of the present invention is that an 
NDGPS system is provided that is more robust, being 
able to detect and recover from equipment failure in one 
of the reference stations. 

65 Another advantage of the present invention is that an 
NDGPS system is provided that compares well to the 
cost of more reference stations and a network of data 
links. 
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These and other objects and advantages of the pres- 
ent invention will no doubt become obvious to those of 
ordinary skill in the art after having read the following 
detailed description of the preferred embodiments 
which are illustrated in the various drawing figures. 

IN THE DRAWINGS 

FIG. 1 is a system block diagram of the worldwide 
differential GPS (WWDGPS) network, and incorporat- 
ing a first embodiment of the present invention; 

FIG. 2 is a block diagram of a typical reference sta- 
tion of the WWDGPS of FIG. 1; 

FIG. 3 is a block diagram of central processing facil- 
ity of the WWDGPS of FIG. 1; 

FIG. 4 is a second level block diagram of the refer- 15 
ence station of FIG. 2; 

FIG. 5 is a data flow diagram for the reference station 
of FIGS. 2 and 4; 

FIG. 6 is a data flow diagram for the central process- 
ing facility of FIG. 3; 

FIG. 7 is a diagram of the orbit position estimator of 
FIG. 6; and 

FIG. 8 is a diagram of the special message Type 25 
transmitted to the users of FIG. 1. 

DETAILED DESCRIPTION OF THE 
PREFERRED EMBODIMENTS 

FIG. 1 illustrates a worldwide differential GPS 
(WWDGPS) network 10. WWDGPS 10 comprises a 
plurality of reference stations 12 in communication with 30 
a plurality of GPS satellites 14 transmitting signals 15 
and belonging to the GPS constellation and a central 
processing facility (CPF) 16 in communication with the 
reference stations 12 b way of a communications satel- 
lite 18 that supports a plurality of very small aperture 
terminal (VSAT) datalinks 19 to CPF 16. In turn, CPF 
16 broadcasts differential GPS correction information 
20 via a communications broadcast satellite 21 to a 
plurality of users 22, which can include DGPS services, 
direct mobile users and even the National Aeronautics 40 
and Space Administration (NASA) in support of the 
well-known space shuttle (STS) missions. A plurality of 
back-up land-links 24 couple the reference stations 12 to 
CPF 16. Candidate sites for the reference stations 12 are 
preferably selected according to the availability of reli- 45 
able power, stable political environments, physical ac- 
cessibility and avoidance of the ionospherically unstable 
regions near the North and South poles. A final selec- 
tion criteria is to prefer international airports, military 
bases and observatories or universities. 

FIG. 2 illustrates a typical reference station 12 which 
comprises a primary GPS receiver/processor 30, a 
backup GPS receiver/processor 32 for redundancy, a 
failure detection and isolation unit 34, a datalink inter- 
face 36, a VSAT satellite link antenna 38, a first backup 55 
datalink 40, a backup VSAT satellite link antenna 42, a 
second backup datalink 44 connected to the landlines 24 
and a meteorological sensor system 48. Receiver/- 
processors 30 and 32 are Trimble Navigation Limited 
model full P-code receivers. Since a significant number 60 
of reference station 12 failures can be caused by site- 
wide catastrophes (e.g., fire or acts of nature), an alter- 
native to including backup elements 32, 40, 42 and 44 is 
to configure WWDGPS 10 with excess reference sta- 
tions 12 and with computational programs in CPF 16 65 
that can tolerate the loss of any one or two reference 
stations 12. Receiver/processors 30 and 32 are sixteen- 
channel L1/L2 P-code survey types. Sensor 12 pro- 
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vides local weather information, such as temperature, 
relative humidity and barometric pressure, to allow 
receiver/processors 30 and 32 to correct for tropo- 
spheric errors associated with the site of the corre- 
5 sponding reference station 12. Reference stations 12 
measure SV data such as pseudorange (PR), accumu- 
lated delta range (ADR), and meteorological data. Tro- 
pospheric and ionospheric delays are computed. Cor- 
rected PR, ADR and ionospheric delays are transmitted 
10 to CPF 16 via VSAT datalink 19 or landlines 24. Data- 
link 19 typically supports 600 bits per second (bps) sim- 
plex data and 150 bps duplex control. 

FIG. 3 illustrates that CPF 16 is comprised of a 
VSAT antenna 50, a primary datalink 52, a backup 
VSAT antenna 54, a first backup datalink 56, a second 
backup datalink 58 connected to the landlines 24, a unit 
60 that processes measurements and corrects and com- 
putes ephemerides, a unit 62 to do system status and 
control functions, a receiver query controller 64, a data- 
20 link processor and formatter 66 and a GPS correction 
broadcast datalink interface 68 to output GPS correc- 
tion information 20. CPF 16 further comprises a central 
processor 69 that is used for general control and execu- 
tion of the items discussed in connection with FIG. 6, 
25 herein. CPF 16 establishes SV integrity by filtering to 
remove reference station clock errors, estimates SV 
clock errors and estimates SV orbit errors. A network 
ionospheric filter produces estimates for worldwide 
ionospheric model corrections. The list of items trans- 
mitted to users 22 comprises an SV integrity message, a 
NAV message validation, SV clock corrections, SV 
orbit corrections and ionospheric model corrections. 

FIG. 4 illustrates a variation on that of FIG. 2 and 
includes measurement pre-processing, in the form of a 
35 process to tropospheric corrections 70 and a process to 
do ionospheric delay corrections 72 that are done 
within a microcomputer 74. Computer-implemented 
processes running on microcomputer 74 can be used to 
implement unit 34 and datalink interfaces 36, 40 and 44. 
Tropospheric delays are useful only to a local user and 
so are not saved. 

FIG. 5 schematically illustrates the data flows in 
reference station 12. Pseudoranges for carriers LI 
(PRL^i) and L2 (PRl 2), ADR and SV clock (Ck) are 
measured by a receiver measurement processor 80. 
PRL1 is output to a tropospheric delay model 82. Ta- 
bles 2A-2C list the input, method and output used to 
implement model 82. Temperature (T), pressure (P), 
and relative humidity (RH) are output from sensors 48 
50 to model 82. Receiver measurement processor 80 out- 
puts PRL1, PRL2 and ADR to an ionospheric delay 
algorithm 84, which produces an ionospheric error 
(E/ono). Tables 3A-3C list the input, method and output 
used to implement model 84. A tropospheric error 
(E tropo) is output to a mixer 86 and is used with E io„ 0 to 
correct a pseudorange error (PRE), and PR-ADR. A 
failure detection and isolation algorithm 88 accepts the 
PRE and PR-ADR to control unit 34. Algorithm 88 
checks the health of each SV 14 and can exclude mea- 
surements that are likely to be erroneous. A communi- 
cations software 90 transmits PRE, Eiono and PR- 
ADR to the communications hardware comprising 
elements 36, 40 and 44. .Tropospheric delay model 82 is 
comprised of one of the following models known to 
those skilled in the art: Chao, Hopfield Zenith, Goad & 
Goodman, Saastamoinen Zenith, Saastamoinen Stan- 
dard, or Davis. Which one to chose is an engineering 
design choice that will be influenced by the user’s needs 
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and the respective advantages and disadvantages of 
each model. 


TABLE 2A 


TROPOSPHERIC DELAY MODEL INPUTS 


£ ij 

— 

elevation angle from station/to S Vj 

K 

= 

orthometric height of station (km) 


*= 

station latitude (degrees) 

P * 

* 

atmospheric pressure (millibars) 

To 

= 

temperature (*K) 

RH 

* 

relative humidity (fraction of 1) 


TABLE 2B 


TROPOSPHERIC DELAY MODEL ALGORITHM 

g' * 1 + 0.0026 cos 2* + 0.00028 ho 

standard gravity 
correction 

15 

[17.157o— 4684] 

water vapor 


e o -6.10*RHe" 

pressure 

20 

Dtz— a00 f7( p,[ 1255 + 0.05 ]e 0 ' 

tropospheric 
^ zenith delay 


c 0.002644 f 0 .UT7 7h n ) 

surface 

25 

u g' 

meteorology 
mapping function 


. G 



X ** Djz 




mapping function 

30 

MF ” . . Ix/n + xll 

Sln + [sin tjj + 0.015] 



t&TROFO,ij = &TZ MF 

code equivalent 
total tropospheric 
delay 

- 35 


TABLE 2C 

TROPOSPHERIC DELAY MODEL OUTPUTS 

AP TROPOJj = code equivalent tropospheric delay for 40 

station, to SVy, at one epoch (time period) 


16 

TABLE 3C-continucd 

IONOSPHERIC GROUP DELAY OUTPUTS 

nation,- & SV}- at one epoch 

CPF 16 receives the reference station 12 measure- 
ments via network communication satellite links. Net- 
work corrections are transmitted to a user 22, which 
may do additional processing, prior to transmitting the 
corrections to another user 22 . A typical user 22 will 
have its own GPS receiver and may have additional 
navigation sensors such as accelerometers, especially 
when used in high dynamic environments. If user 22 has 
an authorized P-code GPS receiver, the SV’s selective 
availability clock dither is removed from the unclassi- 
fied SV clock corrections provided by the CPF 16 . 
Separately, the CPF 16 broadcasts the corrections to 
individual civilian users 22 such as differential GPS 
service providers. Assuming each reference station 12 is 
tracking eight SV’s, each reference station 12 needs a 
communications bandwidth of 300 bps for measurement 
data, 300 bps for data overhead and 300 bps for com- 
mand and control. CPF 16 requires a communications 
bandwidth of 19,800 bps for thirty^three reference sta- 
tions 12 at 600 bps each, NAV data requires 1,200 bps 
for twenty-four SV’s at fifty bps per SV and 9,900 bps 
for command and control 

FIG. 6 illustrates the data flow in CPF 16 . PRE, 
E tono and PR-ADR are received by communications 
hardware 50 , 52 , 54 , 56 and 58 and transferred by a 
communications software 100 to a preprocessor 102, 
which uses a carrier-smoothed, ionospheric-free code 
Hatch filter: 

(/l\ ~fu)] 

Pi, Lb 


Tables 4A-4C list the input, method and output used in 
the implementation of preprocessor 102. 

TABLE 4A 


TABLE 3A 


IONOSPHERIC GROUP DELAY INPUTS 


fl.f2 

= 

Li & L 2 carrier 
frequencies 

Pl1,/> ?L2,ij 


Li & L 2 pseudorange at 
one epoch for i'th 
reference sUtion & S V/ 

Up aij 


elevation & azimuth angle 
for these pseudoranges 
for reference station, & 

_iXi 


PREPROCESSOR INPUTS 



P*,Ll/2 

= 

pseudorange on LI or L2 
(k'th epoch) 

45 

DLL1/2 

— 

Doppler count on LI or L2 
(k'th epoch) 


R,*, R j 

= 

reference receiver (i'th) & 
SV (j'th) position vectors 


fLl» fz.2 

— 

LI and L2 carrier 
frequencies (constants) 

50 

N 

— 

smoothing window maximum 
epoch count 


TABLE 4B 


TABLE 3B 

IONOSPHERIC GROUP DELAY ALGORITHM 


AP JONOjj = 


[P L2,ij - PlI,//] fj,2 

[fii - 


code equivalent ionospheric delay at one epoch 


TABLE 3C 

IONOSPHERIC GROUP DELAY OUTPUTS 

&ij = elevation & azimuth angle for reference 

station, * SV^at one epoch 

AP iONO,ij — line-of-sight ionospheric delay for 


PREPROCESSOR ALGORITHM 
55 



If k > N or there is a cycle slip on the LI or 
60 L2 carrier, set k = 0 and go to the Po^Ll/2 

calculations 

Djt = [f 2 Li/(f 2 Li - f 2 zi>3DjUl ~ *°no- 

[f 2 Z.2/(f 2 X.l - f 2 Z.2)]L>LL2 spheric- 

free 

carrier 

65 m LL 1/2 - 2 D* - D*,JLl/2 code- 

equiv- 

alent, 

integrated 

Doppler 
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TABLE 4B-continued 


5,323,322 


18 


PREPROCESSOR ALGORITHM 


ro.Ll/2 = 


1 


k 

<k+ 1) lJo 


( PlJ.l/2 - ^ 

L 

2 , M m ,z.I/2 

v m=l j 


PiUl/2 - Poxi/2 + J 2 Mjxi/2 


differ- 
ences on 
LI or L2 

carrier- 

smoothed 

initial 

code on LI 
or L2 

carrier- 
smoothed 
code on LI 
or L2 




Apy = P* - [R, - R j) T eij 


'Ll 


Carrier- 

* Smoothed, 

[ (fll - f 2 L2> J KL spheric- 

Free Code 


10 


15 


-HATCH 20 

FILTER 

iono- 

spheric- 

free 

pseudo- 25 

range 

error 




TABLE 4C 



PREPROCESSOR OUTPUT 

k 

= 

carrier-smoothed, ionospheric-free pseudorange 
(epoch*) 

*Pij 


ionospheric-free pseudorange error (Ref. 
Station/, SV/) 


30 


35 


In FIG. 6, ionospheric delays are removed from (PR- 
ADR) and PRE in the preprocessor 80 and the cor- 
rected data is sent to an orbit estimator 104 that uses a 
least-squares and Kalman filter (LS/KF) and a SV 
clock estimator 106 that uses a least-squares (LS) filter. 40 
Tables 5A-5C list the input, method and output used in 
the implementation of orbit estimator 104. An estimate 
iPEPii) is output to an orbit filter stability checker 108. 
Estimator 106 produces a SV clock estimate (e^) that is 
then applied to an SV clock filter and stability checker 45 
110. Filtered zeph * nd e^are sent to a user communica- 
tions interface software 112 and then to a broadcast 
communications software 114 Communications soft- 
ware 100 also provides ionospheric delay measurements 
to an ionospheric estimator 116 with an LS filter to 50 
produce an ionospheric estimate (e ionospheric )• An iono- 
spheric filter stability checker 118 processes the esti- 
mate and sends it to the broadcast communications 
software 114. The estimates are then transmitted out by 
a communication hardware system 120, for example 55 
through interface 68 (FIG. 3). A user communications 
interface hardware transmits out PRE, geph and 
(PR-ADR). 

Orbit estimator 104 is shown in greater detail in FIG. 

7. 60 

The conventional orbit determination approach is to 
use Newton’s equations of motion. An inertial coordi- 
nate system is selected with the origin at the nominal 
center of the earth. Models for all the forces acting on 
an SV 14 are used to obtain the accelerations experi- 65 
enced by S V 14. These are then integrated to obtain the 
orbit of each SV 14. The magnitude of these accelera- 
tions and the corresponding position errors is such that 


not only must the earth, moon and sun gravitational 
force fields be considered for determining SV 14 orbits, 
but also the gravitational forces from ocean and solid 
earth tides, as well as the direct pressure caused by the 
sun and the Earth’s albedo. 

Determining an orbit using Newton’s equations of 
motion that will be fresh enough to last an hour requires 
that the gravitational accelerations from a four-by-four 
earth gravity model, a point mass moon gravity model 
and a point mass solar gravity model be used. For more 
precise orbit estimations, variations in the exact rotation 
and location of the center of rotation of the earth also 
need to be considered, together with the higher order 
accelerations discussed above. The advantage of good 
orbit estimation is that the orbit corrections can be com- 
municated at a substantially lower rate than, for exam- 
ple, the SV clock corrections. The reduced data rates 
are used to lower datalink 19 and 20 requirements. 

The orbit estimator 104 is preferably based on a two- 
step procedure that estimates the orthogonal earth-cen- 
ter fixed orbit position errors for all SV’s at one epoch. 
SV 14 position errors, derived from several epochs, are 
used to estimate the Keplerian element uncertainties of 
the broadcast ephemeris for each Sv. 

Using carrier-smoothed pseudoranges that have been 
tropospherically and ionospherically corrected as in- 
puts, a measurement vector, Y*, can be constructed that 
contains all the pseudoranges at one epoch from all 
reference stations 12. A state vector X* includes refer- 
ence station 12 clock errors, S V 14 clock errors and S V 
14 position errors. A Housholder transformation matrix 
is selected, [Do/?], in which the product of the matrix 
and the measurement matrix [B] of the partials of 
pseudorange errors with respect to the reference re- 
ceiver clock errors, is zero: [Dor][B]= 0. The transfor- 
mation is used to eliminate the unwanted reference 
receiver clock errors in Yfc leading to a corrected mea- 
surement vector y*. A Least-squares estimator is then 
used to compute SV 14 position errors x*for all SV’s at 
one epoch from y*. 

The next step is to aggregate a number of SV position 
error vectors p * from several epoch’s for a single SV, 
into a new measurement vector z. A new state vector 
X 2 is formed that includes the 15 uncertainties in the 
broadcast ephemeris. A least-squares estimator is used 
to estimate these ephemeris uncertainties for the SV 
from z. The procedure is then repeated for all the other 
SV’s. Based on a potential problem in the observability 
of the argument of perigee correction 3o) 0 because to its 
correlation with the mean anomaly SM^it may be desir- 
able to eliminate that error term from the state vector, 
leaving only fourteen SV 14 ephemeris corrections to 
be estimated. 

TABLE 5A 


ORBIT POSITION ERROR ESTIMATOR INPUTS 

Ap/, 

= 

ionospheric-free, carrier-smoothed, pseudorange 

error 



S* 


(SVj visible at ref. station/) 
iine-of-sight unit vector 
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Set up clock and orbit position error observation 
matrices, [B]* and [H]*, based on the SV’s tracked by 
each reference station: 

If the first reference station saw (7) PRN’s: 1 , 4 , ...» 
5 21 ; the second saw ( 5 ): 3 , 4 , . . . , 21 ; and the 33rd saw ( 8 ): 
6 , . . . , 24 ; then: 


[B)k = 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 . 1 . 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

i 

0 

0 

0 

6 

6 . ! . 

0 

6 

0 

0 

6 

_ 0 

0 

0 

0 

0 

0 

o ! 

0 

0 

0 

0 

i 

1 

2 

3 

4 

5 

6 

7 

Station’s 

29 

30 

31 

32 

33 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 . 

1 

2 

3 

4 

5 

6 

7 

PRNs 

20 

21 

22 

23 

24 



o r 

0 r 

0^ 

o r 

or 

0 r ... 

or 

or 

or 

or 

or -i 

0^ 

0^ 

* r i,2 

0 T 

or 

or ... 

or 

or 

or 

or 

or 



or- 

or 

or 

o’r 

o'r ... 

o'r 


o'r 

o’r 

o’r 

0^ 


e lv 

o r 

o r 

or 

or ... 

or 

or 

or 

or 

0^ 

o r 

e r 2,2 

o r 

or 

or ... 

or 

or 

or 

or 

or 


or 

or 

or 

or 

o'r 

o'r ... 

o'r 

4s 

o'r 

o'r 

o’r 


or 

or 

or 

or 

o’r 

o'r ... 

o'r 

o'r 

o'r 

o'r 

* r 33,8 - 

1 

2 

3 

4 

5 

6 

7 

20 

21 

22 

23 

24 


Ap li 

Apn 


[A]k = [B\H] k 

W\fr\y] k ^{D 0 R}k[A\Y] k 

{ &OR }k ~ Combined set 
(vector and matrix) of 
Housholder transforms 


PRN’S 

Composite observation matrix 
Housholder transformed 
observation matrix and vector 


Upy 
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I Ap33^ | 


Check if: [B']*<<[I]? 

If not, write error message or abort orbit position 
estimator 104 for this epoch (time period). 
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21 

Check if: [H']* is full rank? 

If not, write error message or abort orbit position 
estimator 104 for this epoch. 


*k = 


Szi 

6 jc 2 


= 1 mJyk 


5z24 


I k 


zN+n c > 2 zN— dAjy. 


These equations reflect approximations that can be 
5 made because of small eccentricities of the orbits of 
SV’s 14. While the general solution is written in terms 
of sines, cosines and time polynomials, Colombo found 
that for the GPS satellite orbits the specific solution: 

10 

dR(t)=aR CO $inot)+bR sm(noi)+CRt cos {nc^+dRt 
sin (M+e* 


TABLE 5 C 


07tO=^rcos(/i o r)+6rsin(n o /)+C7f cosin^)+dii 
sm(npf) + ej+/7t +^7^ 


ORBIT POSITION ERROR ESTIMATOR OUTPUTS 

x* = orbit position error vector 

(x,y,z; 24 SV's; k'tb epoch) 


FIG. 7 illustrates a method of implementing orbit 
position estimator 104. The method comprises a step 
130 for initializing an epoch counter, a step 132 to do 
orbit position estimation, a step 134 that tests to see if 
the current epoch is the last one, a step 136 that incre- 25 
ments the epoch counter and a step 140 that implements 
an orbit relaxation estimator 142. Preprocessor 102 
feeds in the measurements from SVs 14. Step 132 esti- 
mates the XYZ position for all of the SVs 14 epoch by 
epoch. Each estimate for each SV 14 involves three 30 
parameters. Six to ten epochs are aggregated by step 
134 (in a loop). Each thirty to sixty minutes the aggre- 
gate is supplied to the orbit relaxation estimator step 
140, which estimates fifteen orbit parameter corrections 
each thirty to sixty minutes. The fifteen orbit parameter 35 
corrections are fairly stable and will stay fresh for thirty 
to sixty minutes, so they need only be communicated to 
users 22 that often. 


B N(t) — OR cos (riot) + b& sm(/!o/) + erf cos (n ^) + drf 
sin {n^+eN+frf 

provides good accuracy. Since the solution requires no 
integrations, the unknown coefficients are determined 
by a direct least-squares fit of the GPS measurements 
rotated into orbit plane coordinates. The present inven- 
tion uses least-squares estimating to reduce the compu- 
tational load at various key points. 

Since both the Lagrange Planetary equations and 
Hill’ equations are orbit error equations relative to a 
nominal (broadcast ephemeris) orbit, these approaches 
implicitly assume that the orbit errors are linear. As a 
result these approaches tend to be valid only over short 
orbital arcs and are called orbital relaxation, as opposed 
to orbital estimation, approaches. Tables 6A-6C detail 
the orbit relaxation estimator 142 inputs, method and 
outputs. 

TABLE 6A 

ORBIT RELAXATION ESTIMATOR INPUTS 


CPF SOFTWARE REQUIREMENTS 

A review of the orbit determination literature indi- 
cates that conventional approaches can be separated 
into three categories to formulate the orbit estimation 
problem in terms of the Keplerian (n on-orthogonal) 
elements, the orthogonal orbit plane coordinates, or the 45 
classical inertial coordinates. The differential equations 
describing the orbit errors in Keplerian element form 
are the Lagrange Planetary Equations. The Keplerian 
element formulation is easy to use, because for a per- 
fectly circular orbit, the Keplerian elements are con- 50 
stants. Thus, for the nearly circular orbits of GPS satel- 
lites 14, the changes in the Keplerian elements is small. 
This was one of the considerations by the CPS Joint 
Program Office that led to the selection of the Broad- 
cast Ephemeris in terms of a modified set of Lagrange 55 
Planetary Equations. Specifically, a broadcast ephem- 
eris, which is a best fit of the orbital measurements, is 
made by the system 10 with Lagrange Planetary Equa- 
tions using six initial Keplerian elements, three initial 
Keplerian element rates and six perturbation coeffici- 60 
ents to the Keplerian elements. 

An alternative approach is to use the orbit plane 
orthogonal coordinate approach, the differential equa- 
tions of which are known as Hill’s equations: 

65 

zR-3n 0 2 zR-2n 0 ZT^zAr 


zT+2n 0 Z R = dAr 


Sx*. 8y*, 6z * 

- 

orbit position error 
estimate for one SV 
(k'th epoch) 

r 

— 

radius 

u 

— 

corrected argument of 
latitude 

ft 

= 

corrected longitude of 
ascending node 

I 

== 

corrected inclination 

t k 

= 

kth time point (epoch) 

A 

= 

semi-major axis 

e 

= 

eccentricity 

Ejt 

= 

eccentric anomaly 


- 

argument of latitude 

C U c 

= 

argument of latitude 
correction 

Cn> Crs 

= 

radius correction 

Cjc Cis 

= 

inclination correction 


The second through last items in Table 6A are the 
broadcast ephemeris. 

TABLE 6B 

ORBIT RELAXATION ESTIMATOR ALGORITHM 


Sx* 


PJt = 


£y k 


orbit position error vector for one SV* epoch 


8z* 
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TABLE 6B-continued 
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ORBIT RELAXATION ESTIMATOR ALGORITHM 


Pi 

P2 


general orbit position error vector for one SV 


P k 


5 


10 



*11 

“*12 

“ *13 

*14 

[Tilt = 

*21 

-*22 

*23 

*24 


*31 

“*32 

0 

*34 


where, 

*11 = cos u cos n — sin u cos I sin 11 

ai2 = r [sin u cos H + cos u cos I sin H] 

*13 — r [cos u sin 11 + sin u cos I cos 11] 

*14 — sin u sin I sin 11 

*2i = cos u sin H -j- sin u cos I cos 11 

*22 = r [sin u sin H - cos u cos I cos Hj 

»23 = r [cos u cos n — sin u cos I sin H] 

*24 = r sin u sin I cos 11 

a31 * sin u sin I 

*32 = r cos u sin I 

*34 = r sin u cos I 


15 


20 


25 


30 


[T 2 ]* = 


Bi 

0 

0 

0 


1B12B13B14B15O 0 B18B19O 

0 0 0 B25B26B27O 0 0 

000000000 
0 0 0 B45O 0 0 0 B 4( 10 


0 0 0 0 0 

0 0 0 0 0 

0 0 0 B3.14B3.15 

B 4 ,ll B4.12B4.i3O 0 


35 


where, 


40 


B\\ = —A cos Ek 


B\ g--B 2 7 - 5 4 , 10= cos 24 > k 


A\' 

*2 


M- 


orbit position error sensitivity matrix 


Ah’ 


TABLE 6 B (Continued) 

ORBIT RELAXATION ESTIMATOR ALGORITHM 


*2 = 


8 e 0 

8 A* 

6 n<j 

BM 0 
B <i > 0 

BCus 

BC uc 

BCrc 

BC n 
BCi C 
BC ^ 
Bio 
Bio 
Sfl 0 
8 fi 0 


= {[a] 7 '!*)} -1 W r z 

Orbit relaxation estimator 142 


e* - e + Bc 0 
A' = A + 6A 0 
An' = An + 8n 0 
M' 0 = M 0 + BM 0 
a'o = + Boi 0 


C' uc = C uc -f BC uc 
C'rc — Cn- -f BC k 
C'n = Crs -|- SC ra 
C'/c = C/c -f 8C/c 
C = C^j + BCis 
I o ~ L + 8Ip 
fV<? — + 8H 0 

r = i + Bi 0 
iv = n + sn 0 


Corrected Broadcast 
Ephemeris Parameters 


Bn= 0 -ecosEk) 



TABLE 6C 

B 19- ^^26=^4, 11 = sin 24>* 


ORBIT RELAXATION ESTIMATOR OUTPUTS 


e', 5e 0 


corrected eccentricity & 


— 



correction 

B\3~A e tk/{\-€ cos Ek) 

50 

A', dAo 


corrected semi -major axis & 

B2 5 — (1 + 2 Cyj cos 2 $>k — 2 Cyc sin 24>&) 

An', Bno 

— 

correction 

corrected mean motion 



M o, BM 0 

=3 

adjustment & correction 
corrected initial mean 


Bj 4—A e sin £*/(l -e cos Ek) 
B^ 14= Bi t 12=1 

55 

(ti'o, Bo ) 0 

= 

anomaly & correction 
corrected initial argument 
of perigee & correction 


C us, C'uc BCus, SCus 

= 

corrected argument of 
latitude adjustments & 


£l5 = 2 (C n cos 2 <f>k—C K sin 2 4>*) 


C'rc C'n, BC^ BC n 

= 

corrections 
corrected radius 

Bi, 15 -£4,13-** 

60 

C ic C is, BCjc BCis 

_ 

adjustments & corrections 
corrected inclination 





adjustments & corrections 

•8^4,5— 2 (Cjs cos 2 <j>k-C ic sin 2 4>*) 
and, 


Vo, r, 8Io, Bio 


corrected initial 
inclination, inclination 
rate, & corrections 
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ii'o, n\ an* sh 0 

— 

corrected initial longitude 
of right ascension, 


mi*- mi* [r 2 ] k 




longitude rate of right 
ascension & corrections 
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A network of thirty-three reference stations 12 dis- 
tributed worldwide will satisfy a minimum requirement 
that each satellite 14 can be seen by at least five refer- 
ence stations 12. This minimum number provides the 
raw measurements needed to obtain a high confidence 
position solution for each reference station 12, and aver- 
aged corrections improve with a wider base of measure- 
ments over which to spread the average. In addition, 
the density of the network is sufficient to tolerate a loss 
of two reference stations 12 and still meet a minimum 
visibility requirement of 99%. 

The network of reference stations 12 transmits the 
ionospheric and tropospheric delay corrected GPS 
measurements to the CPF 16. In addition, the measured 
ionospheric delays are also transmitted to the CPF 16. 

An important functional requirement for the network 
is that the density of reference stations 12 be sufficient 
such that each of the twenty-four operational SV’s 14 
can always be seen by at least five reference stations 12. 
This requirement stems from a need to eliminate un- 
wanted reference receiver clock error, in addition to 
estimating SV clock error and three SV position error 
states. 

Each reference station 12 should be housed in an 
environmental enclosure to protect it from precipitation 
and to control the internal temperature from any exter- 
nal temperature extremes, if reference station 12 is not 
located in a manned facility. The GPS antennas must be 
placed outside of the enclosure on a pre-surveyed site 
with unobstructed visibility to the GPS satellite constel- 
lation. Care should also be taken in the placement of 
these antennas to assure that there is a minimum of 
electromagnetic interference and a minimum of multi- 
path reflection from nearby metallic surfaces or bodies 
of water. 

TROPOSPHERIC DELAY CORRECTION 

The tropospheric delay in the GPS measurements is 
an unwanted delay in the pseudorange and accumulated 
delta range data that is introduced by the troposphere. 
The delay, which is RF frequency-independent, typi- 
cally introduces a delay of two meters for straight 
above head vertical (zenith) line-of-sight elevation an- 
gles. The delay increases to about twenty-five meters 
for line-of-sight elevation angles of 5°. The total delay 
consists of a dry and a wet component. The dry compo- 
nent accounts for about 90% of the total delay and the 
wet component accounts for the rest. While the dry 
component can be reasonably well modeled without 
any meteorological sensor data, the wet component 
requires measurements of the local weather conditions 
along the line-of-sight to satellite 14 for maximum accu- 
racy. Typically, surface meteorological sensor measure- 
ments are extrapolated to the higher altitudes, using 
conventional atmospheric profiles. 

Any of several tropospheric delay models can be 
used. The ideal algorithm provides the best worldwide 
accuracy, estimates both the dry and wet components 
and has the largest set of input variables. The algorithms 
with a large input variable set can be expected to pro- 
vide better worldwide all-weather performance, over 
algorithms with a smaller input variable set. A compari- 
son of these and other tropospheric delay models rela- 
tive to ray-trace tropospheric delay estimates was per- 
formed by Janes, H. W., et al, in “Analysis of Tropo- 
spheric Delay Prediction Models: Comparisons with 
Ray-Tracing and Implications for GPS Relative Posi- 
tioning (A Summary)”, GPS'90, Ottawa, Canada, Sep- 
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tember 1990. It was determined that the Saastemoinen 
zenith delay model combined with the Davis mapping 
function and the Goad and Goodman water vapor algo- 
rithm is a model to be preferred. The principal inputs to 
5 the tropospheric delay model are temperature, pressure 
and relative humidity measurements from the meteoro- 
logical sensor module 48. The relative humidity is then 
converted to water vapor pressure, which is used to- 
gether with the temperature and pressure measurements 
10 to determine the tropospheric zenith delay. Line-of- 
sight elevation angles, are also needed from GPS re- 
ceivers 30 and 32. In addition the geodetic latitude and 
altitude of the GPS antenna is required, which are es- 
15 tablished at the time the antenna for GPS receivers 30 
and 32 is surveyed at its permanent location. This data 
is then used to determine the tropospheric delay map- 
ping function as well as the standard gravity correction 
in the zenith delay model. 

20 IONOSPHERIC DELAY CORRECTION 

The ionized components of the atmosphere, which 
extend from about fifty kilometers high to 1,000 kilome- 
ters, will introduce delays in RF signals that vary in- 
25 versely as the square of the RF frequency. At GPS LI 
and L2 frequencies, the delays can approach fifty me- 
ters for low elevation angle lines-of-sight. The delay 
slows down the pseudorange (code) measurement and 
advances the accumulated delta range (carrier) mea- 
30 surement. The ionosphere has a variable total electron 
count (TEC) that determines the extend of the delay. 
The TEC is variable with respect to local time and the 
geomagnetic latitude. In addition to daily (diurnal) vari- 
ations, ionospheric delay varies with the month of the 
35' year and the eleven year solar sunspot cycle. Iono- 
spheric delay is practically unpredictable in the polar 
regions, around the equator and everywhere during 
solar storms. 

P-code reference receivers are needed for reference 
40 stations 12 because they have the ability to use dual-fre- 
quency measurements to fmd the ionospheric delay 
directly, using a simple algorithm that only depends on 
the carrier frequencies of and the time delay between 
the two measurements. 

45 For C/A code receivers, the Navigation Message 
(NAV) from each GPS satellite contains a model for the 
ionospheric delay that was developed by Klobuchar. 
The Klobuchar model is a fit to the more extensive Bent 
model, which is based on a comprehensive data base of 
worldwide ionospheric measurements. The Klobuchar 
model fit is best in the northern hemisphere and in par- 
ticular, in the mid-latitudes over the continental United 
States. It is formulated as a truncated four parameter 
55 cosine model and includes a night time bias, an ampli- 
tude (four coefficient polynomial), a constant phase and 
a period (four coefficient polynomial). These ten coeffi- 
cients, which are included in the Navigation Message, 
are based on monthly averages. It is estimated that a 
60 broadcast ionospheric delay model corrects for about 
50% to 60% of the total ionospheric delay. 

The inputs to the ionospheric group delay algorithm, 
above, are pseudorange measurements from the LI and 
L2 channels, in addition to the accumulated delta range. 
65 The algorithm computes the ionospheric delay along 
the line-of-sight to each of the SV’s in view. The iono- 
spheric delays are removed from the pseudorange error 
estimates and are also sent to CPF 16, separately. 
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FAILURE DETECTION AND ISOLATION 
ALGORITHM 

Both receivers 30 and 32 operate at the same time 
with one receiver designated as a lead, or primary, 5 
receiver. The measured data from both receivers 30 and 
32, as well as additional data, such as signal strength, is 
fed to failure detection and isolation algorithm 88. The 
algorithm will check the integrity of the data being 
collected by the lead receiver. If the integrity of the 10 
data from the receiver is in question, the data is flagged 
to CPF 16, which may then elect to switch to the 
backup receiver 32, using a command link from CPF 16 
to reference station 12. 

REFERENCE RECEIVER 15 

In an experiment conducted by the inventors, a dual- 
frequency, unauthorized P-code survey receiver was 
selected. Originally a twenty-channel receiver was se- 
lected as the nominal receiver, since it provides ten 20 
channels of tracking on each frequency. A Monte Carlo 
simulation was done to examine the number of opera- 
tional SV’s visible to a reference receiver located any- 
where in the world. For a 10“ mask angle, no more that 
eight SV’s are visible 99% of the time. Thus, a sixteen- 25 
channel receiver is adequate and can be used as the 
preferred reference station receiver. In building an 
error budget for reference station 12, some measure- 
ments will exhibit an accuracy that varies between 0.6 
meters for zenith line-of-sight measurements and 1.5 30 
meters for line-of-sight measurements at a 10° elevation 
angle. A decision to use unauthorized P-code receivers 
should be based on understanding that P-code is en- 
crypted (Y-code) only intermittently. 

A full Y-code GPS receiver is available to classified 35 
military users only. An unauthorized full P-code re- 
ceiver has a full-wavelength codeless tracking mode as 
a backup. An example of that type of receiver is the 
Rogue SNR-8 receiver that reverts to a codeless cross- 
correlation tracking mode when Y-code is turned on. 40 
Other full P-code receivers can revert to codeless track- 
ing mode based on squaring the L2 channel, for exam- 
ple, the Ashtech P-12 receiver. An L2 P-code receiver 
reverts to a codeless tracking mode based on squaring 
the L2 channel. The Trimble Geodetic Surveyor IIP is 45 
an example of that type of receiver. The Trimble Geo- 
detic Surveyor II is an example of dual -frequency, re- 
ceiver that lacks P-code and tracks the L2 channel by 
squaring the channel. The agreement in the delays be- 
tween a P-code and a codeless tracking mode is good 50 
for elevation angles higher than 30° and deteriorates for 
the lower elevation angles for the codeless tracking 
mode. When using carrier phase derived delays as a 
reference for the code delays, good agreement exists 
between the code and carrier delays when a receiver is 55 
operating in the P-code tracking mode. Code delays 
deteriorate for the lower elevation angles when the 
receiver is operating in the codeless tracking mode. It 
can be concluded that the absolute L1-L2 delay are 
obtained, whether a receiver is operating in the full 60 
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a noisy reference signal for the unknown P-code signal 
on the second frequency to obtain the relative P-code 
delay, without decoding the unknown P-code signals. 

The choice of what receiver type to select is a com- 
plex one. The receiver type that may be adequate for a 
test reference receiver, which is checked in a non-Y 
Code environment, but it may not be adequate for an 
operational reference receiver. Also, the GPS receiver 
technology is still evolving such that a receiver selected 
as reference receiver may be outdated by the time an 
operational worldwide differential GPS network 10 can 
be implemented. 

METEOROLOGICAL SENSOR 

Meteorological sensor measurements are required to 
estimate a “wet” component of the tropospheric delay, 
which accounts for about 10% of the total tropospheric 
delay. While more sophisticated sensor packages can be 
selected to measure the tropospheric delay directly 
along the line-of-sight to each SV, surface measure- 
ments of temperature, pressure and relative humidity 
are good enough for the operational reference station 
12. The preferred tropospheric delay algorithm extrap- 
olates these surface measurements through the tropo- 
sphere using standard atmospheric profiles. A generic 
meteorological sensor module 48 could be assembled by 
purchasing the individual sensors and the data acquisi- 
tion hardware and software components and integrating 
these into a custom module. The sensor voltages, from 
the temperature, pressure and relative humidity sensors 
are preferably low-pass filtered, to remove any noise. 
These filtered voltages are then digitized to obtain the 
raw pressure, temperature and relative humidity mea- 
surements, in engineering units. Since the pressure and 
relative humidity sensors are temperature-sensitive, 
they must be corrected for any temperature-caused 
errors, unless the sensors are temperature controlled. 
The corrected meteorological measurements are then 
passed to the tropospheric delay correction algorithm 
82. 

An alternative solution to assembling a custom mete- 
orological sensor module 48 is to select a portable 
weather station with sensors of the required accuracies. 
An error budget for the tropospheric delay compensa- 
tion hardware and software should be based on the 
quoted accuracies of the algorithm and a linear error 
analysis for the algorithm to establish the sensor accu- 
racy requirements. The limiting accuracy is due to the 
mapping function at low elevation angles. Thus, sensors 
that are more accurate than the mapping function are 
not needed. Based on this, pressure and temperature 
sensors with an accuracy of 1 % to 2% each and relative 
humidity sensors with an accuracy of 5% to 10% 
should be adequate for use in meteorological sensor 
system 48. A 1% pressure sensor has a 10 millibar (mb) 
uncertainty and a 1% temperature sensor has a 5“ Fahr- 
enheit (F) uncertainty relative to absolute zero. Other 
factors that need to be considered are the environmental 
specifications, the robustness of the design and hard- 
ware and the amount of testing that has been done. 


P-code or the codeless tracking mode. However, in a 
codeless tracking mode, the L1-L2 delay becomes nois- 
ier for the lower elevation angle measurements. This 
stems from the fact that P-code tracking loops use a 
receiver generated noiseless reference signal to decode 65 
the received P-code signals when not encrypted. In the 
codeless cross-correlated tracking mode of the receiver, 
the unknown P-code signal on one frequency is used as 


REFERENCE PROCESSOR 

Reference station microprocessor 74 hosts the tropo- 
spheric and the ionospheric delay algorithms 82 and 84 
and communications software 90. A 286-type IBM-AT 
compatible personal computer or clone can be used for 
microcomputer 74. The ionospheric delay algorithm 84 
and possibly the tropospheric delay algorithm 82 can be 
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installed as firmware in the reference receiver proces- 
sor. A principal function of microprocessor 74, is to 
host failure detection and isolation software 88, in addi- 
tion to any communications software 90. 

5 

NETWORK COMMUNICATIONS SYSTEM 

A key function of the network communications sys- 
tem comprising datalinks 19 and satellites 18 is to trans- 
mit the measurements from the worldwide network of 
reference stations 12 to CPF 16 accurately and with 1 
minimum latency at a minimum of cost. The communi- 
cations system must also permit control of remote, unat- 
tended reference stations 12 by CPF 16. Because of the 
worldwide distribution of reference stations 12, the 
satellite communications datalink 19 is required. The 
International Telecommunications Satellite Organiza- 
tion (INTELSAT) is able to provide the digital data 
service required. The service uses the inexpensive Very 
Small Aperture Terminals (VSAT’s) and is fully inter- 20 
active. Both “Ku” and “C” band transponders with 
global and hemispherical beams are available. Band- 
width allocations of up to 72 MHz are also available. A 
link 19 from each reference station 12 to CPF 16 is 
preferably a full-time, non-preemptive lease from IN- 25 
TELSAT. The use of the INTELSAT service requires 
the use of three INTELSAT communications satellites, 
in order to obtain the necessary global coverage. Since 
CPF 16 will probably not be able to simultaneously 
communicate with all three INTELSAT satellites, 30 
feeder links (or alternately multiple central processing 
facilities 16) may be necessary. In case a primary VS AT 
communications link 19 is lost, a second (backup) satel- 
lite communications link is a good idea. The first choice 35 
for the backup link is a second INTELSAT service 
using different satellites. The same VSAT’s could be 
used for both the primary and the secondary link, in that 
case. The secondary link is not necessarily required on 
a full-time basis. An alternative backup satellite link 40 
could comprise the International Maritime Satellite 
Organization (INMARSAT) Standard “A” or Standard 
“C” service. For reference stations 12 located in areas 
that offer a public switched telephone system, such a 
system can be used as a backup land-line communica- 45 
tions system, as in landlines 24. 

A number of communications link performance issues 
need to be considered. Two key issues are the minimum 
latency of the received data and the data integrity. 
While there is some transmission delay of network data 50 
over the satellite links 19, the use of a global system of 
reference stations 12 needing to communicate with a 
single CPF 16 imposes a requirement to use a satellite 
communication system. If the delay is a critical issue, ^ 
multiple CPF’s 16 can be used, instead. Data integrity in 
transmitting the data over the communications links 19 
are assured by redundant communications link interface 
hardware at reference stations 12 and CPF 16 (e.g., 
through the use of redundant VSAT’s 18). It can also be ^ 
improved by reliable protocols and an insensitivity of 
the communications signals to weather. 

Since usage cost is a very important operational fac- 
tor, tradeoffs need to be made between data volume 
versus usage costs. Data volume is reduced by selecting 65 
only key data for transmission, by using data compres- 
sion techniques (such as truncating the data accuracy to 
that actually required), and by using mixed data rates. 


REFERENCE STATION COMMUNICATIONS 
LINK SOFTWARE 

A key function of reference station 12 communica- 
tions link software 90 is to format reference station 12 
data to be sent to CPF 16 in messages for transmission 
over the communications link 19 Formatting includes a 
proper scaling or truncation of the data from reference 
station 12, in addition to a scheduling of data that is 
transmitted at different data rates. The software 90 must 
also decompress or decode control messages received 
from CPF 16. The signal-to-noise ratio is included to 
check the health of the SV 14, as well as to monitor any 
multipath signal. If any data is bad and outside its nomi- 
nal dynamic range, it is sent to CPF 16 after having 
been scaled. 

The Navigation Message includes the broadcast 
ephemeris, while the issue of data ephemeris (IODE), 
identifies the age of that ephemeris. Only one Naviga- 
tion Message is required from each SV114 from the 
whole network. Since each SV 14 is tracked by seven 
reference stations 12, there would otherwise be a seven- 
fold duplication of Navigation Messages, if each refer- 
ence station 12 sends the Navigation Message from all 
S V’s in view. Thus, to assure that only one set is sent for 
each SV, CPF 16 could request specific reference sta- 
tions 12 to send the Navigation Message from specific 
SV’s 14. Alternately, a protocol can be set up so that 
each reference station 12 sends back the Navigation 
Message for any SV 14 that has a local elevation angle 
above some minimum value. 

DATA TRANSMISSION REQUIREMENTS 

Data transmission requirements dictate the required 
communications link 19 bandwidth. The bandwidth 
estimate for the data to be sent between reference sta- 
tion 12 and CPF 16 are derived from the number of bits 
of data that need to be sent, multiplied by the required 
sampling frequency for the data. It is possible to esti- 
mate the communications bandwidth required to trans- 
mit the data from reference station 12 to CPF 16. For 
eight SV’s 14 being tracked at each reference station 12, 
approximately 600 bps is required, including a 300 bps 
overhead to allow for headers, parity bits, etc. The total 
does not include the Navigation Message, at fifty bps, 
for those eight SV’s. Only twenty-four Navigation Mes- 
sages need to be transmitted from a network of thirty- 
three reference stations 12, which is less than one Navi- 
gation Message, at fifty bps, per reference station 12. 
Control messages are sent back to reference station 12 
by CPF 16. Approximately twenty-one kbps of data is 
received by CPF 16 from the thirty-three reference 
stations 12. The total includes the Navigation Messages 
for the 24 SV’s. About 9.9 kbps of command and con- 
trol bandwidth are sent by CPF 16 back to the thirty- 
three reference stations 12. 

CPF COMMUNICATIONS LINK SOFTWARE 

The principal function of CPF 16 communications 
link 19 software 90 is to receive and decode the data 
sent from the thirty-three reference stations 12 and to 
format command and control messages back to the 
individual remote reference stations 12. Due to the 
volume of data, a key issue is how to schedule the ar- 
rival of the data from the Network without data colli- 
sions or loss of data. While this is not an issue for the 
prototype CPF 16, it is an issue that requires additional 
analysis for the operational system. 
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In FIG. 8, a special message 150 is constructed by 
broadcast software 114 to economize in the volume of 
data placed on datalinks 20. CPF 16 produces a Radio 
Technical SC-104 format Type 25 message to supple- 
ment pseudorange corrections (PRC) transmitted in 5 
Type 1 messages. In system 10, the PRC’s in the Type 
1 message apply to a user 22 located at coordinates that 
are transmitted in a Type 3 message to CPF 16 by user 
22 over datalink 20. The correction applicable to a user 
22 is: 10 

PR C* PR Cf\ + {bt u — latTti'pLat-ni + ( ion u — 
bn n>pLon n5 

where, PRC T\ is the PRC from the Type 1 message in meters, latu Wld 
lon u are the latitude and longitude of the user 22 in 13 
radians, lat/3 and Ion 73 are the position from Type 3, 
pLatj25 and pLonns are the latitude and longitude 
correction factor from Type 25 in meters per radian. 
FIG. 8 shows the configuration of special message 150, 
which comprises a thirty-bit word 152 containing 24 20 
bits for a header word one and a six bit parity; a thirty- 
bit word 154 containing 24 bits for a header word two 
and a six bit parity. A block for a first satellite 16 com- 
prises a thirty-bit word 156 containing three bits of fill, 
five bits of satellite identification, sixteen bits of Type 25 25 
latitude correction factor in meters per radian; a thirty- 
bit word 158 containing six bits for issue of data, sixteen 
bits of Type 25 longitude correction factor in meters per 
radian. A variable number of blocks each containing a 
pair of words 160 and 162 follow and are similar to 30 
words 156 and 158, respectively. An “n /A ” block con- 
taining a pair of words 164 and 166 comprises a thirty- 
bit word 156 containing three bits of fill, five bits of 
“n^” satellite identification, sixteen bits of Type 25 
latitude correction factor in meters per radian; a thirty- 35 
bit word 158 containing six bits for issue of data, sixteen 
bits of Type 25 longitude correction factor in meters per 
radian. 

CENTRAL PROCESSING FACILITY ^ 

A primary function of CPF 16 is to take the iono- 
spherically and tropospherically corrected measure- 
ments from the reference stations 12 and estimate the 
SV clock and ephemeris errors. The unknown reference 
receiver clock errors are not corrected for, but rather 45 
are considered nuisance variables that are removed in 
the preprocessing stage prior to the estimation of these 
SV errors. (The variable is removed by estimating the 
SV clock.) Another function of CPF 16 is to provide a 
filtered estimate of the ionospheric delay for users 22 50 
that do not have dual-frequency GPS receivers. The 
ionospheric estimate is based on an ionospheric filter or 
estimation algorithm that uses the dual-frequency GPS 
reference receiver ionospheric measurements transmit- 
ted to CPF 16 from reference stations 12. 55 

Since reference stations 12 are at remote locations 
and are usually to be unattended, another function of 
CPF 16 is to query the status of reference stations 12. If 
needed, CPF 16 will then send a required command to 
the remote receivers. 60 

EXTERNAL INTERFACES 

Principal inputs to CPF 16 are reference station 12 
measurements that have been received over the com- 
munication system. As required, CPF 16 will also assert 65 
control over the remote unattended reference stations 
12 via the communication system. As shown in FIGS. 2 
and 3, the communication system interface uses a 
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VS AT both for a primary and the backup interface. A 
land line interface is also provided as another backup 
interface, where the interface are implemented. 

CPF 16 interfaces with a user 22 over datalink 20. 
The specific form of the interface depends on the spe- 
cific user. However, as discussed above, the interface 
also includes the option for VSATs and land line com- 
munications interfaces. Users 22 receive the SV clock 
and ephemeris corrections, the ionospheric delay model 
estimates and SV integrity information. 

SV ORBIT ESTIMATOR STABILITY CHECKS 

Estimates obtained from orbit estimator 104 are 
checked for stability of the filter solution. This includes 
checks done on the least-squares residuals to determine 
that there are no significant biases in these residuals. 
Other checks should also be done, as required. 

SV CLOCK ESTIMATOR 

SV clock estimator 106 is presented in FIG. 6. Tables 
7A-7C show the inputs, method and outputs of SV 
clock estimator 106. The algorithm uses carrier- 
smoothed pseudorange errors, which have been cor- 
rected for tropospheric delays, ionospheric delays and 
SV orbit errors. These adjusted pseudorange errors 
from all the reference stations 12 at the same epoch are 
stacked together to form a measurement vector, Z. A 
state vector b is formulated that includes the reference 
receiver clock errors on top and S V 14 clock errors on 
the bottom. To remove unwanted reference receiver 
clock errors, a difference operator [Dce] is used. Use of 
the matrix on the measurement matrix [F] which con- 
sists of the partial derivatives of the pseudorange errors 
with respect to the reference receiver clock errors, 
eliminates the matrix [Dce][F]= 0. Following the oper- 
ation, a least-squares estimator is used to estimate S V 14 
clock errors, 6T, for all the SV’s, using the corrected 
measurement vector z. 

TABLE 7A 


SV CLOCK ESTIMATOR INPUTS 


Pij^pk 


ionospheric-free pseudorange 

R / 

= 

i'th reference station position vector 


= 

line-of-sight unit vector (station,-, SVy) 

Ml! 


measurement noise (diagonal) variance 

e' 

= 

corrected eccentricity 

A' 

35= 

corrected semi-major axis 

An' 

as 

corrected mean motion adjustment 

U'o 

s= 

corrected initial mean anomaly 

a>'o 

= 

corrected initial argument of perigee 

C'hJ, C'j lc 

* 

corrected argument of latitude 
adjustments 

C'm 

= 

corrected radius adjustments 

C'ic C'ir 


corrected inclination adjustments 

\o 

= 

corrected initial inclination 

I' 

= 

corrected inclination rate adjustment 

ft'o 


corrected initial longitude of right 
ascension 

O' 

“ 

corrected longitude rate of right 
ascension 

P 

* 

WGS84 value of earth’s universal 
gravitational parameter 



GPS system time at transmission & 
ephemeris reference epoch 

O f 

= 

WGS84 value of earth’s rotation rate 


TABLE 7B 

SV CLOCK ESTIMATOR ALGORITHMS 


Corrected Ephemeris Calculations 
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TABLE 7B-continued 


SV CLOCK ESTIMATOR ALGORITHMS 


p 

. + An' 

(A') 3 


toe 

M'* - M'o + t* n' 

E'* = M'* + e' sin E'* 

Solve iteratively for E'* 
with initial estimate for 
E'* as M'*, on right side 
of equation. 


f * = arctan 


NTT7T 


sin E'k 


cos E'* — e' 


4>'k = f A + W'o 

Au k — C x/s sin 2(f> k +- 0'^ c 

cos 2<f>’k 

Ar'* = C' re cos 24>'* + C'„ 
sin 2<p'k 

AI'* - C’i c cos 2<f>'* + C'is 
sin 2 <j>k 

u'* - 4>* + Au'* 

r'* = A' [1 - e' cos E'*] + Ar'* 
I*Jt — L 0 + AT* 4- t* I’ 
ft'* = ft' 0 + t* ft'o “ t ft e 

x' k = r'* [cos u'* cos ft'* - 
sin u'* cos I'* sin ft'*] 
y’ k = r'* [cos u'* sin ft’* + 
sin u'* cos I'* cos ft'*] 
z’* = r'* sin u* sin I'* 


R '/ = 


corrected mean motion 


time from ephemeris 
reference epoch 
corrected mean anomaly 
corrected eccentric 
anomaly 


corrected true anomaly 


corrected argument of 
latitude 

corrected argument of 
latitude adjustment 

corrected radius 
adjustment 

corrected inclination 
adjustment 

corrected argument of 
latitude 

corrected radius 
corrected inclination 
corrected longitude of 
ascending node 
corrected 


SV position 
(earth-fixed) 


10 


15 


20 


25 


-10000 
0 0 0 -1 0 


-continued 

o o 
o o 


ooooo 

0 0 0 0 0 


0000000... 0 -1 000 

00 -1 0000. ..ooooo 

000 -1 000... ooooo 


0 0 0 0 0 0 0 


0-1000 


0000000 . 0000-1. 
1 2 3 4 5 6 7 ... 20 21 22 23 24 

PRN^s 


z = 


Pce] = 


30 


35 


40 


45 


50 


1 


Ap'n 

Ap '12 


A P'ij 
Ap'33j 

-1 


measurement vector at epoch r* 


N 2Af 1 1 N 2 M 22 

1 1 

M 63/11 \ 6M22 M 123/33 

1 1 1 


0 

-2 


0 

-3 


MlUJn N i2M 22 M 123/33 N 123/44 
1 


N/(/-/)A/ll 


0 

-ALzlL 


U 


87 = 


S7i 

87 2 


8724 


Z = [Z>ce] z 
[G’] = [Dce) [Grj 


- {[G] T [G'}}-' [G') T z 

SV Clock Estimator 

TABLE 7C 


SV CLOCK ESTIMATOR OUTPUTS 


SV CLOCK ESTIMATION CALCULATIONS 


55 


t* 

8T 

R'; 


epoch 

SV clock corrections at the epoch 
corrected SVy position at the epoch 


A p'y = pij— [R R r j\ T tij corrected pseudorange error 
Set up SV clock sensitivity [G] matrix based on SV’s 
tracked by each reference station, 
e.g.: If the first reference station saw PRN’s: 1, 4, . . . 
, 21; the second saw 3, 4, . . . 21; and the 33rd saw: 
6, . . . , 24; then the [G] matrix looks as follows: 

[G] = 


SV CLOCK ESTIMATOR STABILITY CHECKS 

The results of SV clock estimator 106 are checked to 
assure a stable solution. This will probably include a 
check of the least-squares estimator residuals to check 
for whiteness or lack of any significant biases. 

IONOSPHERIC MODEL ESTIMATOR 

The function of the ionospheric model estimator 116 
is to provide ionospheric delay estimates for any civilian 
user who is using a single-frequency GPS receiver and 


60 


65 
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wishes to obtain better accuracy than obtained with the 
broadcast ionosphere model, the Klobuchar Model. 
Tables 8A-8C show the inputs, methods and output of 
ionospheric model estimator 116. 

In its simplest form, the ionospheric delay can be 5 
characterized in terms of a zenith delay and an obliquity 
factor, or mapping function, which translates the zenith 
delay to non-vertical lines-of-sight from the user to SV 
14 being tracked. 

Klobuchar formulated a model of the ionospheric 10 
delay which uses a zenith delay that is a truncated co- 
sine function. It also uses a mapping function that de- 
pends only on the line-of-sight elevation angle and as- 
sumes that the ionosphere is concentrated in a shell at a 
fixed altitude of 350 kilometers. The four principal pa- 15 
rameters used to describe the Klobuchar model are: a 
night time bias (DC) and a day time amplitude that is 
multiplied by a cosine function that has a period and a 
phase. Klobuchar selected the order of the polynomial 
for these four parameters. Thus, he selected a night time 20 
bias (O’th order polynomial) that might have been more 
accurately represented with a fourth order polynomial. 
Klobuchar approximated the day time amplitude with a 
fourth order polynomial. Klobuchar selected a constant 


36 

comparable accuracy and is reputed to be less computa- 
tionally complex, the model and its revision (IRI-90) are 
good candidates for use as ionospheric model estimator 
116. The preferred approach is based on deriving a 
contour VTEC map, using the dual-frequency iono- 
spheric delay measurements at reference stations 12. In 
its simplest level, the dual-frequency ionospheric mea- 
surements are rotated to the equivalent local (sub-ionos- 
pheric shell penetration of the line-of-sight) zenith de- 
lay. These zenith delays are then fit to a two-dimen- 
sional surface for the whole globe. A mapping function, 
as well as regional values of that zenith delay surface 
are then incorporated by user 22 to determine its line-of- 
sight ionospheric delays. The specific ionospheric delay 
model equation for the preferred approach is: 

D\{rXhj) *= 




phase of 1,400 hours, although a fourth order polyno- 
mial is more accurate. Klobuchar selected a fourth 
order polynomial for the period, although a higher 
order polynomial can be more accurate. Klobuchar 
represented the cosine function as a four-term series 
expansion, rather than using the cosine function itself or 
considering higher order terms. The reason for making 
these approximations were dictated, in part, by the 
desire to keep the broadcast ionosphere model rela- 
tively simple for the single-frequency user as well as to 
minimize the amount of data to be sent over the GPS 
navigation message. 

A summary of some of the more sophisticated global 
ionospheric models indicates whether the electron pro- 
file and the integrated electron profile, or TEC, is avail- 
able in these models. It also shows the input variables 
that the models accept; hopefully the more variables 
that are input, the more accurate the solution that is 
obtained. The ITS-78, IONCAP and Bradley models 
cannot be used for the ionospheric model estimator 116 
since these models do not extend completely through 
the ionosphere. Of the remaining models, the Bent and 
IRI-86, International reference Ionosphere of 1986, are 
the preferred models. The recently released update to 
the IRI-86 model, the IRI-90 model, may also be con- 
sidered. 50 

There are a number of options open for enhancing the 
broadcast ionospheric delay model. The ten parameters 
used in the Klobuchar model can be estimated directly. 
Better accuracy can probably be achieved since hourly 
measurements can be used rather than monthly aver- 55 
ages to determine these parameters. Another option 
consists of estimating perturbations, in the form of a bias 
and a scale factor, to the broadcast model. Finally, 
another option is to use higher order expansions for the 
four parameters of the Klobuchar model as well as to 60 
evaluate the cosine function directly, rather than using 
the truncated Taylor series expansion. While all of these 
approaches have their merits, they are all limited by the 
fact that the cosine model used is inherently not accu- 
rate enough. 65 

Another approach is to use the Bent model in various 
ways. The Bent model appears to be too complex to be 
used as a real-time model. Since the IRI-86 model is of 


The empirical model uses either a two-dimensional 
polynomial expansion, or alternately, a two-dimensional 
Fourier series in terms of geomagnetic latitude and local 
time. The mapping function is also approximated using 
a first order Taylor series about the nominal ionospheric 
altitude. The variability of the ionospheric altitude can 
be expressed either as a two-dimensional polynomial 
expansion or as a two-dimensional Fourier series in 
geomagnetic latitude and local time. Dual-frequency 
measurements are used in the equation to derive the 
unknown coefficients of the zenith delay and iono- 
spheric altitude variation. 

A closer look at the candidate ionospheric delay 
model indicates that the estimation problem is still non- 
linear. This arises from the fact that the zenith delay 
expansion and altitude variation expansions occur as a 
product. The unknown coefficients, however, can be 
estimated iteratively. Specifically, the zenith coeffici- 
ents are estimated first, assuming the altitude variation is 
negligible. Next, using the initial estimate of the zenith 
delay, the altitude effect is estimated. The results of the 
first pass are then compared to the raw measurements 
and the residuals are checked to see if another iteration 
is required. 

While the model is driven directly by the dual-fre- 
quency measurements, a reference model, such as the 
IRI-90 model can be employed during a startup phase to 
select either the polynomial or the Fourier series expan- 
sions as well as to determine the order of these expan- 
sions. 

TABLE 8A 

GLOBAL IONOSPHERIC DELAY MODEL INPUTS 


AP JONOij 


measured line-of-sight 
ionospheric delay (station/, 
SV,) 

X|, $,• 

= 

latitude and longitude of 
i'th reference station 

a ij 

- 

station / - SV,- line-of-sight 
elevation and azimuth angle 

T GPS 

= 

GPS time 

hio 

— 

nominal (thin shell model) 
ionospheric altitude 

Xat, 

= 

latitude and longitude of 
geomagnetic north pole 
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TABLE 8A-continued 

GLOBAL IONOSPHERIC DELAY MODEL INPUTS 

*£ — spherical earth radius 

[, J(I) — number of stations & number 

of SV’s observed at station/ 

M = . zenith ionospheric delay 

expansion limits with respect 
to normalized local time and 
local geomagnetic latitude 

U K. — ionospheric altitude 

variation expansion limits 
with respect to normalized 
local time and local 
geomagnetic latitude 

-mean. Cmu = ionospheric model mean and 

rms residual convergence test 
limits 


TABLE 8B 

GLOBAL IONOSPHERIC DELAY MODEL ALGORITHM 
Calculate Independent Variables 
Ah nj — 0 initial estimate of ionospheric altitude 
variation 

h nj — h jo 4- Ah/y local (subionic) ionospheric altitude 
A ij = (90 0 - aj) - sin- 1 {(R£Cos c,y)/(R£ + h/,y)} central 
earth angle: site — subionic point (deg) 

X/y = sin -1 {sin X/ cos Ay + 

cos X/ sin Ay cos a,y} local 

(subionic) geocentric latitude 

Bjij — Bf 4- sin - 1 {sin Ay sin ay/cos X/y} local 

(subionic) longitude 

X ; y = sin - 1 {sin X/,y sin A 3/ 4- 

cos X/y cos Xm cos (Biij - 

0m)} local (subionic) geomagnetic latitude 

ty — 240 Bjij 4- Tqps local time 


240 

(ty - 86,400) 
240 

(t,y 4- 86,400) 
240 


for normalized local time. 

I AP/CWOl.l I 

AP/CW 01.2 


0 < r,y < 86,400 


t/» > 86,400 
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-continued 

i/.Xr,X)= k la hlk T'\ k 

1=0 k= 0 

previous (polynomial) estimate of ionospheric altitude 
variation from nominal, or 


AAXt.X) = flAOO + ^2 i [anok cos(2 kk) 4- ahoj(+k sin(2 k\)] + 


L a h 10 cos(1t) 4- 

2 4- 

/-I flA£+l,0 sin(lr) 


cos(1t)cos(2A:X) 4- 

L K ahL+ijc+kCOs(\Tysm(2k\) + 
/=1 *=1 0A 2 £ + i, 2£ + * sin(lr)cos(2AX) 4- 

CA3£+1,3£+* sin(lr)sin(2^X) 


previous (Fourier series) estimate of ionospheric alti- 
tude variation from nominal. 1 
2 5 For initial pass, Ah/=0 


J _ f (/?£ COS Cy) 2 

L (^£ + */G>) 2 _ 


off-nominal ionospheric altitude scaling coefficient 


tfllCn 

0 

0 . 

0 

0 

0 

CnBn 

0 . 

0 

0 

0 

0 

0 . 

CuBij 

, 0 





. 0 

0 

0 

0 . 

0 

ClJQpIM 


scaling matrix 

(Matrix [S] 2 , below, is too wide to fit across the page, 
and so is split into two sections . . . ) 


Dl ~ kPjONOiJ 


line-of-sight measured ionospheric 
delay vector 


X11 

X 2 .i 

■ x»„ 

X12 

X 2 .2 • . 

• • 

X U{I) 

X 2 V(/) 



r io ? c %% ER 


IONOSPHERIC ZENITH DELAY ESTIMATOR 
(POLYNOMIAL/FOURIER SERIES) 


• nominal mapping function 


Til • T^ll TllXn . T^nX^n 

T 12 - T N n T12X12 • ^12X^12 

T U(l) ■ ^IrKD T V+(Mw 


zenith (polynomial) delay ionospheric 
measurement sensitivity matrix 


(Re 4- h\o) 2 


(The following matrix is so wide it is split into six 
sections.) 


[53 z = 
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-continued 


1 

cos( 2 Xn) 

cos ( 4 Xn) 

cos (2JfXii) 

sin( 2 Xn) 

sin ( 4 Xn) 

I 

cos ( 2 Xn) 

cos (4Xl2) 

cos ( 23 /Xl 2 ) 

sin (2X12) 

«n (2Xi2) 

1 

cos (2X47(7)) 

cos (4X747)) . . 

. cos 

«n (2X/44) 

tin (4X747)) 


sin( 2 AfXn) cos(th) cos( 2 ru) 

sin cos (ru) cos (2rn) 


cos (Nru) nn(rn) tin (2th) 

cos (Nr 12) tin (n 2 ) sin (2ti 2 ) 


sin ( 23 /X/^t)) cos (r/^/)) cos (2 r/^/)) 


cos sin (r/^) sin (2r/^/)> 


... sin (Nr n) 

. . . sin (Nr 12) 


cos (r n ) cos ( 2 Xn) 
cos (ri 2 ) cos (2X12) 


cos (rn) cos ( 4 Xn) 
cos (T12) cos (4X12) 


cos (Nr n) cos (2JAn) 
COS (Nri2) cos (2MX12) 


... sin (Ntj^d) cos cos (2X7^7)) 


cos (r/yg)) cos (4X/^/)) 


cos (Nt/^;)) cos (lAfki^f)) 


cos (th) sin ( 2 Xn) 
cos (ri 2 ) sin ( 2 Xj 2 ) 


cos (rn) sin ( 4 Xji) 
cos (ri 2 ) sin ( 4 Xi 2 ) 


cos (Nr 11) sin ( 2 ATX 11) sin (rn) cos ( 2 Xu) 

cos (Nr 12) sin ( 2 M\n) sin (ri 2 ) cos ( 2 Xi 2 ) 


cos sin (2X7,47)) cos (r/^/)) sin (4X447)) 


cos (Ntjad) sin (2MX447)) sin (T447)) cos (2X747)) 


sin (r n ) cos ( 4 X U ) 
sin (rj 2 ) cos ( 4 X ]2 ) 

sin (r/^/)) cos (4A747)) 


sin (Nr 11) cos ( 2 MXn) sin (th) sin ( 2 Xn) 

sin (Nri2) cos (IWX12) sin (rj 2 ) sin ( 2 Xn) 


sin (Nr /^T)) cos (23/X747)) sin ^7475) sin (2X47(7)) 


sin (r u ) sin ( 4 Xn) 
sin (rj 2 ) sin ( 4 Xj 2 ) 

sin (T747)) sin (4X747)) 


sin (Nrn) sin ( 2 A/Xn) 
sin (Nri 2 ) sin ( 2 AfXj 2 ) 

sin (Nr 7,7(7)) sin (2JWX747)) 


35 


coo 

«01 


C 0 , 2 M 

c l,0 


40 


a 2 = 


fl 2N,0 

oil 

flj 2 




The above is a zenith (Fourier series) delay iono- 
spheric measurement sensitivity matrix. 

The zenith delay measurement matrix is: 

[#]=[C][S] 


a nm 


a 4N,AM 


The ionospheric (polynomial) zenith delay coefficient 
estimator is either a: 


50 


2 ? 7 Z(T f X) = 2 2 a nm T n \ m ionospheric zenith 

(polynomial) delay 
estimate for r and X 


coo 

«01 


<*z 


a nm 




a NM 


55 M 

= ooo+ 2 km cos( 2 roX) + 00*7+ m sin(2mX)] + 
m= 1 


N 

X [o„o cos(nr) + c* + ,,,o sin(rtr)] + 
n=l 


60 


N M 

X X [o„ m cos(nr)cos( 2 mX)] 
n= 1 m=l 


m cos(nr)sin( 2 mX) + « 2 N+ji, 2 J#+m »n(nr)cos( 2 mX) + 

65 fl 37 vr +n , 3 j»/ +m sin(Hr)sin( 2 mX;], 


or an ionospheric (Fourier series) zenith delay coeffici- The above is the ionospheric zenith (Fourier series) 
ent estimator: delay estimate for r and X. 
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-continued 


SDjj - { Dizirij, \ij) Bjj - A PjoNOij) line-of-sight 
ionospheric delay residual 


/ = 2 J(f) total number of measurements collected 
i=l 

by network at one epoch 

SDmean = *4“ 2 P BD„ residual mean 
J i= 1 j—l 


Aij — (1 — Qy) Bij 


A\\ 0 0 

0 A\2 o 


M- 


10 


0 0 
0 0 


0 

0 

0 

0 


scaling matrix 


= 


1 / J(f) , 

. tt 2 2 (SDjj) 2 root-mean-square residual 
A J i=l j— 1 

If {(8D mean < Cmean) and (SDrms < Crms r)}» 

convergence test passed (exit ionospheric delay model). 

IONOSPHERIC ALTITUDE VARIATION 
ESTIMATOR 


6£>n 

6^12 


15 


ionospheric (polynomial) altitude sensitivity matrix: 


BDm 


bd l 


8D im 


[S\h = 

1 x„ . , 

x 2 n 

1 x, 2 . . 

X 2 12 

■ ■ .. 
Til 

..X' 12 .. 

T12 

• • T^nTiiXn . 

• T^12 
T 12^12 

■ . rL n * K n 

■ . t L 12^. K U 


1 ^u(n ■■ 

x IAI) 

■ x k 1A d . . 

^IAI) 

■ ■ 

or 

(Matrix [S]/r is so wide it is split into six sections.) 


ionospheric delay measurement 
residual vector 




1 COS ( 2 Xn) 

1 cos (2X12) 

cos ( 4 Xn) 
cos ( 4 X j2 ) 

cos ( 2 AXn) 
cos ( 2 AXi 2 ) 

sin ( 2 Xn) sin ( 4 Xu) 

sin (2X12) sin (4X12) 



1 cos (2 X/,y ( /)) 

cos ( 4 X/ f /( j)) 

... cos (2AX/ f y(7)) sin ( 2 X/^j)) sin ( 4 X/,/(j)) 



sin ( 2 A"Xj j) 
sin ( 2 AX 12) 

cos (r n) cos( 2 r U ) 

cos (rj 2 ) cos ( 2 ti 2 ) 

cos (Lrn) 
cos (Lr 12) 

sin (tu) 
sin (ti 2 ) 

sin ( 2 th) 
sin ( 2 ti 2 ) 


sin ( 2 fCkj >J{T )) 

cos cos (2t/j ( /)) 

... cos sin (t/^(/)) 

sin { 2 t jjw) 


sin ( Lt u ) 
sin (Ltu) 

cos(ru) cos ( 2 Xn) 
COS (TJ2) cos (2X12) 

cos (tu) cos ( 4 Xn) 
cos (tj 2 ) cos ( 4 Xi 2 ) 

. . . cos (Ltu) cos ( 2 AXn) 
... cos (Lr 12) cos (2AX12) 


sin (Lr/^/)) 

cos (JU(T)) cos (2X^(7)) 

cos (t/^/)) cos (4 

• • • cos (Ltjj^j)) cos ( 2 AX/ r /(y ) ) 

cos(tii) sin ( 2 Xn) 
cos (ti 2 ) sin (2X12) 

cos (Til) sin ( 4 Xn) 
cos (T12) sin ( 4 Xi 2 ) 

. . . cos (Ltu) sin ( 2 AXn) 
cos (Lt 12) sin ( 2 AXi 2 ) 

sin (tji) cos (2Xji) 
sin (T12) cos (2X12) 


cos (TIJV)) sin (2 X/v(i)> cos (r/^j)) sin (4 \jjy)) cos sin {IKKj^) sin (r/^) cos 


sin (rn) cos (4X U ) 
sin (7*12) cos (4X12) 


sin (tj M ) cos (4X/ t y { i)) . 


sin {Lt\ j) cos (2/CXn) 
sin (Lr 12 ) cos (2 AXj 2 ) 


sin (Lt/jw) cos 


sin (rn) sin (2Xn) 
sin (rn) sin ( 2 X 12 ) 


sin {TiJiD) sin (2 


sin (rn) sin(4X n ) 
sin (t 12 ) sin (4 Xi 2 ) 


sin (tj M ) sin (4 Xjj(i)) 


sin (Lrn) sin (2 AXn) 
sin (Lt 12 ) sin (2AXi 2 ) 


sin (Ltj^d) sin 
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The above is an ionospheric (Fourier series) altitude 
sensitivity matrix. 


[Hlb=[A)[S\ h 

ionospheric altitude measurement matrix, 


5 


0* 00 
0*01 


Oh 


= {mh) T mh)}-Him T w 


OHLK 


10 


15 


ionospheric (polynomial) altitude variation coefficient 
estimator 20 


0*00 

0*01 


25 


0*0, IK 
0*10 


0 * - 


0A2Z..O 

0*11 

0*12 




30 


0 * 1 ,* 


35 


0*41,4* 


ionospheric (Fourier series) altitude variation coeffici- 
ent estimator. Return to ionospheric zenith delay esti- 
mator. 


TABLE 8C 


GLOBAL IONOSPHERIC DELAY MODEL OUTPUTS 



= 

zenith ionospheric delay model 
coefficient vector 

*A 

= 

ionospheric altitude variation 
model coefficient vector 

SD meart> 

- 

ionospheric model mean and 

SDmu 


rms residuals 


The ionospheric model estimator 116 is driven by 
dual-frequency GPS reference receiver (30 or 32) mea- 55 
surements of the observed absolute ionospheric delays. 
The estimator is a global model, although a user algo- 
rithm may be based on a local subset, or truncated ver- 
sion, of the global model. Because of the time varying 
nature of the ionosphere (the reasons for which are not 
well understood), the model is updated hourly. Al- 
though it is desirable to have the model be as accurate 
as possible, the global extent of the model with a reason- 
able number of coefficients and the updating these coef- 
ficients only once an hour will probably limit the accu- 
racy to 80-90% of the ionospheric delay. Also, the need 
to minimize communication costs and realize a compu- 
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tationally simple user algorithm will further limit the 
accuracy at the user end. 

IONOSPHERIC FILTER STABILITY CHECKS 

The specific form of the filter stability checks re- 
quired for ionospheric model estimator 116 will depend 
on the estimator. The stability checks are based on 
checking the estimator residuals for whiteness or any 
significant biases. 

SYSTEM STATUS AND CONTROL FUNCTION 
SOFTWARE 

Since reference stations 12 are in remote locations 
relative to CPF 16, they are unmanned, and it is desir- 
able to have the capability to command the two refer- 
ence receivers in each reference station 12. The data 
from the two reference receivers is compared by failure 
detection and isolation software 88 at reference station 
12. With one of the receivers (30 and 32) designated as 
primary and the second as backup, the designations can 
be switched if it is determined that the data from the 
primary receiver is suspect. The quality control can be 
done by failure detection and isolation software 88. The 
authority to switch from the primary to the backup 
receiver 32 is handled automatically at reference station 
12, with a manual over-ride capability from CPF 16. 
Alternately, the authority is fully assigned to reference 
station 12 or to CPF 16. 

CENTRAL PROCESSOR 

The choice of an operational central processor 16 
requirement is strongly influenced by the code require- 
ments of SV 14 orbit estimator 104, SV clock estimator 
106 and the ionospheric model estimator 116. Most 
likely a work station environment is required, since a 
PC environment will not be adequate. 

For an operational central processor 69, a number of 
different options can be considered. The output display 
requirements for the operational central processor 69 
are dictated by the need for operator control and moni- 
toring of the various software modules as well as the 
capability to provide WWDGPS network situation and 
status displays. These dual requirements could be satis- 
fied with a single large screen monitor with windowing 
capability to display both network status as well as 
central processing status. Alternately, an operator con- 
sole consisting of a standard monitor and a separate 
large screen display to present the network status could 
be employed. 

BROADCAST COMMUNICATIONS SYSTEM 
FUNCTIONS 

The principal function of datalink 20 is to transmit S V 
14 clock and ephemeris corrections and the network 
ionospheric model corrections to users 22 and can in- 
clude companies that sell differential corrections to 
their customers or end-users of these corrections. 

The corrections must be broadcast in near real-time 
and some of the corrections need to be sent out more 
frequently, to minimize data latency. A two-sigma error 
growth occurs because of selective availability as a 
function of the delay in getting the corrections to user 
22, assuming the corrections are perfect. For a thirty 
second delay, nearly five meters (2cr) error may be 
introduced by SA clock dither. In order to minimize the 
cost of using datalink 20, the corrections must be pack- 
aged efficiently. 
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BROADCAST COMMUNICATIONS SYSTEM 
DATA FLOW REQUIREMENTS 

The section discusses the data flow requirements for 
datalink 20 by presenting the message definition and the 
transmission requirements. The broadcast message is 
separated into two integrity messages, a primary integ- 
rity data and a navigation message validation. It also 
includes three correction messages: a ephemeris correc- 
tion data, a clock correction and a ionospheric correc- 
tion data. 

The primary integrity data specifies the integrity of 
all the satellites. The navigation message validation, in 
turn, specifies the integrity of the navigation message. 
The ephemeris correction data includes in addition to 
the ephemeris corrections, estimated by SV 14 orbit 
estimator 104 , the time of the message, the issue of data 
ephemeris (IODE), and SV ID. The clock correction 
data includes in addition to SV 14 clock corrections, the 
time of the message, SV 14 id and SV 14 signal quality 
or signal strength, because of its importance, the pri- 
mary integrity data is broadcast at one Hz. because of 
the SA dither on SV 14 clock, the clock corrections 
should be broadcast at a 0.1 Hz frequency to user 22. 

While the ionospheric correction frequencies depend 
in part on the ionospheric model estimator 116 , these 
probably are comparable to the frequencies required for 
correcting SV 14 orbit errors. The civilian message 
transmission requirement is under 160 bps. conserva- 
tively, it is estimated that data rates of 300 bps should be 
adequate for all the data, excluding the ionospheric 
data. The ionospheric data can require an additional 300 
bps for a total of 600 bps. 

BROADCAST COMMUNICATIONS 
HARDWARE REQUIREMENTS 

It is anticipated that datalink 20 will have two levels. 
The first level uses communications satellites and the 
data is broadcast using global or hemispherical beams. 
The level will cover the entire globe from about 70° 
south to 70° north latitude. The second level uses local 
area communications. There are two systems that can 
be used for the first level. 

The first is INMARSAT COMSAT, which is satel- 
lite system that has a data service using small, inexpen- 
sive terminals and omnidirectional antennas. The ser- 
vice is intended for two-way communication but can be 
used in the broadcast mode. It has the capability to 
broadcast data to users 22 at 600 bps. 

The second service, the INTELSAT INTELNET 
service, has two modes. One is a point-to-point service 
and the other is a point-to-multipoint service. As previ- 
ously discussed, the communication system uses the 
point-to-point service. The point-to-multipoint service, 
as normally configured, requires the use of small direc- 
tional antennas. However, changes are made to the 
system so that only an omnidirectional antenna is re- 
quired and allows the system to still meet the 600 bps 
data rate requirement. 

The second level of datalink 20 is exemplified by use 
of a subcarrier from an FM broadcast station. The sub- 
carrier broadcast could transmit the data to just a local 
area and users 22 would require a communication re- 
ceiver (FM radio) much less expensive than that re- 
quired (satellite receiver) by users of the first level of 
datalink 20. 
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REDUCING DATA COMMUNICATIONS 
VOLUMES 

The broadcast ephemeris parameter information 
being sent to users 22 by CPF 16 is preferably done at a 
rate not so frequent that the broadcast communication 
datalink 20 becomes congested, about ten minutes is the 
minimum period. However, the rate must not be so 
infrequent that it is impossible to obtain a reasonable set 
of ephemeris parameters that will retain their accuracy. 
Two hundred minutes is the limit. 

Data communications volume and bandwidth in a 
networked differential GPS system can be reduced by 
deriving ionospheric correction data (PRCs) in relation 
to a master reference position x"/ 1 obtained from a 
plurality of fixed location reference stations. Interpolat- 
ing said PRCs to produce ionospheric correction data 
for a plurality of user positions s uwr amongst said refer- 
ence stations. Communicating partials of said PRCs to 
said users wherein each user can adjust received differ- 
ential corrections for its own location by applying a 
simple linear function using measured gradients as coef- 
ficients to difference between user position x^ er and 
master reference position x 7 */ 1 : 

. DPR C /3 gradient of PRC in northerly direction, 

a PRC/dE~ gradient of PRC in easterly direction, and 

PRC user =^ PR C^f + (dPR C/dN, d PRC/d E, 

0)'(x user —x re fl). 

INTERPOLATED CORRECTIONS FOR USERS 

Three or more reference stations 12 are used to find 
averaged corrections for arbitrary points around the 
world that lie amongst reference stations 12. These 
points produce what appears to be a plurality of phan- 
tom reference stations to the users 22. Corrections that 
are sent to users 22 are normalized (referenced to) to the 
nearest phantom reference position, in two-dimensions, 
East-West and north-south. Unlike the prior art, the 
accuracy of corrections provided to users 22 is not a 
function of individual separation distances of a user 22 
from a particular reference station 12 or a phantom 
reference station. Each reference station 12 contributes 
to a blend of weighted averages that have a consistent 
accuracy at substantially all points on the globe. The 
two-dimensional partials with respect to the phantom 
reference station points are compensated (fudged) to 
ease the computational load on users 22 and result in a 
high quality correction. Users 22 are not uniform in 
their position solution accuracy needs, so those users 
that do not require extreme accuracy are sent shorter 
(less bit precision) messages. Less precise messages are 
less expensive to produce because the computational 
loads can be scaled back in real-time and accordingly 
subscription prices could be reduced. 

RAPID SATELLITE ORBIT EPHEMERIDES 
PARAMETER DETERMINATION 

Above, the orbit ephemeris parameter estimator re- 
lies on a history of SV positions collected over several 
epochs. Using SV carrier measurements that have been 
corrected for ionospheric and tropospheric effects, SV 
velocity can be obtained directly. A rapid GPS network 
SV orbit ephemeris parameter estimator that receives 
both SV position and velocity such that an immediate 
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Kepler is available is comprised of means for obtaining 
measurements from a plurality of reference stations 
comprising carrier-smoothed pseudoranges that have 
been corrected for both ionospheric and tropospheric 
delays and that produces an ionospheric delay measure- 
ment. Such a means is reference station 12 (FIG. 1). A 
central processing facility, such as CPF 16, is used for 
computing SV positions from said reference station 
measurements. A carrier correction means for correct- 
ing an S V carrier for both ionospheric and tropospheric 
delays is needed at each reference station 12. This pro- 
vides a precise measurement of corresponding SV ve- 
locity that is communicated to CPF 16. Included are 
estimating means for generating a set of SV orbit 
ephemeris parameters which processes all of said refer- 
ence station measurements through one of a least- 
squares estimator or Kalman filter, such as has been 
described above for CPF 16. A set of SV orbit ephem- 
eris parameters is then broadcast to users 22. The accu- 
racy of navigation of each user 22 is substantially im- 
proved thereby. 

Although the present invention has been described in 
terms of the presently preferred embodiments, it is to be 
understood that the disclosure is not to be interpreted as 
limiting. Various alterations and modifications will no 
doubt become apparent to those skilled in the art after 
having read the above disclosure. Accordingly, it is 
intended that the appended claims be interpreted as 
covering all alterations and modifications as fall within 
the true spirit and scope of the invention. 

What is claimed is: 

1. An ionospheric delay modeling method, the com- 
prising the steps of: 

obtaining both code and carrier measurements from a 
dual-frequency GPS receiver of ionospheric delays 
in a sparse network of reference station receivers 
tracking SV’s in the GPS constellation; 

using a two-dimensional geomagnetic latitude and 
local time Taylor Series or Fourier Series to map 
equivalent local zenith measurements into an iono- 
spheric delay global model; and 

using a two-step iterated least-squares estimator that 
estimates zenith coefficients for assumed iono- 
sphere altitude, which checks the goodness-of-fit 
and examines residuals, which estimates ionosphere 
altitude for assumed zenith coefficients and that 
repeats the foregoing steps when required to 
achieve convergence. 

2. The method of claim 1, further comprising the 
subsequent step of: 

truncating estimates to obtain a regional model from 
said ionospheric delay global model. 

3. The method of claim 1, further comprising the 
subsequent step of: 

local-fitting to obtain a regional model from said 
ionospheric delay global model. 

4. A differential GPS (PGPS) system for providing 
GPS correction information to a plurality of users ac- 
cording to each user's position, comprising: 

a regional network of a plurality of reference stations 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV with unauthorized 
GPS receivers; 

a central processing facility (CPF) in communication 
with the regional network of reference stations and 
having means for computing satellite orbit and 
clock errors and self-synchronizing clock means 
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for deriving a DGPS system clock from a plurality 
of said unauthorized GPS receivers; and 
communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users, wherein 
the use of complex ionospheric and tropospheric 
delay models are avoided and authorized GPS 
receivers are unnecessary. 

5. The system of claim 4, further comprising: 
ionospheric delays computing means connected to 

the CPF for computing ionospheric delays for 
single-frequency users; and 
ionospheric message broadcasting means connected 
to the CPF for broadcasting to particular users 
separate ionospheric messages. 

6. The system of claim 5, wherein 

said ionospheric messages comprise near real-time 
estimate of parameters tuned to particular regions 
of the earth. 

7. The system of claim 4, wherein: 

each of the plurality of reference stations have at least 
one GPS receiver for estimating ionospheric delays 
which are communicated to the CPF for determi- 
nation of the ionospheric delays in a region of the 
earth by interpolation of a plurality of GPS signals 
received by all of said GPS receivers. 

8. The system of claim 4, further comprising: 
group processing means connected to the CPF for 

reducing communications costs at the reference 
stations having means for separation of the refer- 
ence stations into groups, each group having means 
for independently measuring satellite positions and 
clocks; and 

group control center connected to the CPF for input- 
ting collective estimates for SV orbits such that a 
reduction in communications and data rates are 
realized. 

9. The system of claim 8, wherein: 

each group control center has an atomic clock for 
providing time synchronization between the 
groups of reference stations. 

10. The system of claim 8, wherein: 

each group control center comprises means for qual- 
ity control and robustness improvement based on 
atomic clocks. 

11. The system of claim 4, wherein: 

the plurality of reference stations is limited in number 
by a predetermined minimum number of reference 
stations needed to keep all SV’s in said GPS con- 
stellation of satellites in view by at least one refer- 
ence station at any one time wherein re-initializa- 
tions of any particular SV are avoided. 

12. The system of claim 4, further comprising: 
broadcast communication means connected to the 

CPF for sending to users DGPS corrections with 
differing accuracy levels with an identical message. 

13. The system of claim 4, wherein the CPF further 
comprises: 

measurement collection means for obtaining both 
code and carrier measurements from a dual-fre- 
quency GPS receiver of ionospheric delays over a 
period not longer than the period an individual SV 
is in view from a single reference station in a sparse 
network of reference station receivers tracking 
SV’s in the GPS constellation; 
measurement preprocessing means for preprocessing 
measurements to remove tropospheric and mea- 
sured ionospheric signal delays; 
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ephemeris estimating means for estimating SV orbit 
ephemeris parameters using one of a least-squares 
estimator and Kalman filter; and 
ephemeris broadcasting means for broadcasting SV 
orbit ephemeris information to a plurality of users. 5 

14. The system of claim 13, wherein the means for 
broadcasting ephemeris information to users comprises: 

ephemeris information transmission means for send- 
ing said ephemeris information at a rate not exceed- 
ing once every ten minutes and not less than once 10 
every two hundred minutes. 

15. The system of claim 13, wherein the measurement 
collection means comprises one of: a P-code GPS re- 
ceiver, a cross-correlation codeless GPS receiver, or a 
combination P-code and cross-correlation codeless 15 
GPS receiver. 

16. The system of claim 13, wherein said SV orbit 
ephemeris information sent by the ephemeris broadcast- 
ing means comprises: 

a complete set of SV orbit parameters. 20 

17. The system of claim 13, wherein said SV orbit 
ephemeris information sent by the ephemeris broadcast- 
ing means comprises corrections to an SV orbit parame- 
ter set that has been downlinked by an individual SV. 

18. The system of claim 13, wherein the ephemeris 25 
estimating means comprises: 

orbit position estimating means for estimating a set of 
all SV orbit positions in aggregate at each epoch 
from said preprocessed measurements using a least 
squares estimator, wherein one estimator is used for 30 
all SV’s; and 

orbit parameter estimating means for estimating SV 
orbit parameters for each SV from an accumulated 
history of positions for each SV using at least one 
of a least squares estimator and a Kalman filter, 35 
wherein a single estimator is used for each SV. 

19. The system of claim 18, wherein the means for 
estimating a set of all SV orbit positions comprises: 

clock error removal means for removing SV and 
reference station clock errors using a linear trans- 40 
formation of the preprocessed measurements; and 
SV position estimating means for estimating SV posi- 
tions from the transformed measurements using a 
least squares estimator. 

20. The system of claim 13, further comprising: 45 

means for removing SV and reference site clock er- 
rors using an orthogonal transformation; 

means for estimating SV orbit position errors for all 
SV’s at one epoch using a least-squares estimator; 
means for repeating said preprocessing, removing 50 
and estimating every few minutes until sufficient 
measurement estimates have been obtained; and 
means for propagating all position estimates for one 
SV to a common epoch. 

21. The system of claim 13, wherein said CPF further 55 
comprises: 

means for using carrier-smooth pseudoranges that 
have been tropospherically and ionospherically 
corrected as inputs to construct a measurement 
vector Y*that contains all of a plurality of pseudo- 60 
ranges at one epoch from a plurality of reference 
stations. 

22. A dual-frequency differential GPS network 
(DGPS), comprising: 

a plurality of reference stations not exceeding thirty- 65 
three in number substantially equally distributed 
around the earth at locations on land for continu- 
ous tracking of every one of the satellites in the 


GPS constellation by at least one of the reference 
stations at any one time; 

a central processing facility (CPF) in communication 
with all of the reference stations; 
a plurality of dual-frequency non-authorized GPS 
receivers wherein at least one each is located at 
each of the reference stations and provide measure- 
ments of signals from various GPS satellites for the 
computation of ionospheric delays in the local area 
of the respective reference station; and 
computing means in communication with the plural- 
ity of reference stations for estimating global iono- 
spheric delays for a plurality of user locations and 
for computing satellite orbit and clock corrections 
without authorized precision code access. 

23. The DGPS of claim 22, further comprising: 
communication means connected to the computing 

means for supplying interpolated ionospheric de- 
lays and satellite orbit corrections to particular 
users in relation to the locations of said users. 

24. A method of reducing data communications vol- 
ume and bandwidth in a networked differential GPS 
system comprising a plurality of reference stations, a 
central processing facility (CPF) and group control 
centers, the method comprising the steps of: 

estimating GPS satellite orbit ephemerides in four- 
dimensions and in real-time at a CPF or a group 
control center; 

communicating GPS satellite clock corrections de- 
rived from the step of estimating at a first update 
rate from the CPF or group control centers to a 
plurality of users; and 

communicating GPS satellite orbit corrections de- 
rived from the step of estimating at a second update 
rate from the CPF or group control centers to a 
plurality of users, said second update rate being 
substantially slower than said first update rate. 

25. A method of producing real-time GPS satellite 
corrections in a networked differential GPS system, the 
method comprising the steps of: 

receiving a plurality of GPS position signals with a 
plurality of codeless GPS receivers distributed 
amongst a plurality of reference stations; 
processing said GPS position signals with a GPS 
navigation processor to determine orbit position 
and orbit ephemerides for at least one orbiting GPS 
satellite; 

estimating GPS satellite orbit position errors and 
their orbit ephemerides from said GPS positions 
signals with a least squares estimator; and 
estimating GPS satellite clock errors from said GPS 
position signals with a least-squares estimator. 

26. A method of estimating GPS satellite orbit param- 
eters to produce stable estimates that stay fresh for at 
least ten minutes, the method comprising the steps of: 

receiving a plurality of GPS position signals with a 
GPS receiver; 

processing said GPS position signals with a GPS 
navigation processor to determine orbit position 
and orbit ephemerides for at least one orbiting GPS 
satellite; 

initializing an epoch counter; 
estimating the XYZ orbit position of a plurality of 
GPS satellites received by said GPS receiver at one 
epoch; 

repeating the estimating of the XYZ orbit position for 
at least ten minutes; and 
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feeding said orbit position estimates gathered over an 
epoch to an orbit relaxation estimator which esti- 
mates at least fifteen orbit parameter corrections 
each thirty to sixty minutes. 

27. A rapid GPS network space vehicle (SV) orbit 5 
ephemeris parameter estimator that receives both SV 
position and velocity such that an immediate Kepler is 
available, comprising: 

means for obtaining measurements from a plurality of 
reference stations comprising carrier-smoothed 1° 
pseudoranges that have been corrected for both 
ionospheric and tropospheric delays and that pro- 
duces an ionospheric delay measurement; 
a central processing facility for computing SV posi- 
tions from said reference station measurements; 15 
carrier correction means for correcting an SV carrier 
for both ionospheric and tropospheric delays at 
each reference station thereby providing a precise 
measurement of corresponding SV velocity; 
estimating means for generating a set of SV orbit 20 
ephemeris parameters which process all of said 
reference station measurements through one of a 
least-squares estimator or Kalman filter; and 
means for broadcasting said set of SV orbit ephemeris ^ 
parameters to a plurality of users such that user 1 
navigation is substantially improved thereby. 

28. A differential GPS system, comprising: 

a plurality of reference stations at separate locations 
and each having a codeless dual-frequency GPS 
receiver for obtaining ionospheric and tropo- 
spheric delay measurements at more than one of 
said reference stations from a single GPS satellite; 
a central processing facility in communication with 
the plurality of reference stations and having means 35 
for synthesizing clock error signals from non-preci- 
sion code measurements communicated from said 
codeless dual-frequency GPS receivers. 

29. A GPS network space vehicle (SV) orbit ephem- 
eris estimator, comprising: 4 q 

means for obtaining both code and carrier measure- 
ments from a dual-frequency GPS receiver of iono- 
spheric delays over a period not longer than the 
period an individual SV is in view from a single 
reference station in a sparse network of reference 45 
station receivers tracking SV’s in the GPS constel- 
lation, wherein each reference station comprises 
means for ionospheric processing with a codeless 
L2 receiver to use night-time measurements that 
represent stable initial estimates of ionosphere and 50 
for using ionospheric rates from L2 to get estimates 
of ionospheric delays at zenith; 
means for preprocessing measurements to remove 
tropospheric and measured ionospheric signal de- 
lays; 55 

means for estimating SV orbit ephemeris parameters 
using one of a least-squares estimator and Kalman 
filter; and 

means for broadcasting SV orbit ephemeris informa- 
tion to a plurality of users. 60 

means for estimating SV orbit ephemeris parameters 
using one of a least-squares estimator and Kalman 
filter; and 

means for broadcasting SV orbit ephemeris informa- 
tion to a plurality of users. 65 

30. A differential GPS (DG PS) system that provides 
GPS correction information to a plurality of users ac- 
cording to each user’s position, comprising: 


52 

a global network of a plurality of reference stations 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV; 
a central processing facility (CPF) having means for 
computing satellite orbit and clock errors, compris- 
ing: 

means for using carrier-smoothed pseudoranges 
that have been tropospherically and ionospheri- 
cally corrected as inputs to construct a measure- 
ment vector Y* that contains all of a plurality of 
pseudoranges at one epoch from a plurality of 
reference stations; 

means for selecting a linear transformation matrix, 
[Do/?], in which the product of the matrix and a 
measurement matrix [B] of partials of pseudo- 
range errors with respect to reference receiver 
clock errors, is zero:[Do/?][B]=0, said linear 
transformation being used to eliminate unwanted 
reference receiver clock errors in Y* which leads 
to a corrected measurement vector y*; and 
a least-squares estimator for computing SV posi- 
tion errors x* for all SV’s at one epoch from y*; 
and 

communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users. 

31. A differential GPS (DGPS) system that provides 
GPS correction information to a plurality of users ac- 
cording to each user’s position, comprising: 

a global network of a plurality of reference stations 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV; 
a central processing facility (CPF) having means for 
computing satellite orbit and clock errors; and 
communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users; 
wherein each of the plurality of reference stations 
include means for using measurements from the 
SV’s to estimate a zenith component of ionospheric 
delay at each reference station; and 
the CPF comprises means to combine said zenith 
components to calibrate a model of ionospheric 
delay as a function of time and comprises means to 
distribute said model to a plurality of single-fre- 
quency users. 

32. The system of claim 31, wherein: 

said zenith component estimating means comprises a 
least-squares estimator. 

33. A differential GPS (DGPS) system that provides 
GPS correction information to a plurality of users ac- 
cording to each user's position, comprising: 

a global network of a plurality of reference stations 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV; 
a central processing facility (CPI 7 ) having means for 
computing satellite orbit and clock errors; and 
communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users; 
wherein each reference station comprises means for 
the measurement of slant ranges that are mapped to 
reference station coordinates. 
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34. A differential GPS (DGPS) system that provides 
GPS correction information to a plurality of users ac- 
cording to each user’s position comprising: 
a global network of a plurality of reference stations 5 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV; 
a central processing facility (CPF) having means for 
computing satellite orbit and clock errors; and 10 
communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users; 
wherein the CPF comprises means for generating 15 
quality control information by using over-specified 
reference station information and applying a root- 
sum-of-squares algorithm. 
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35. A differential GPS (DGPS) system that provides 
GPS correction information to a plurality of users ac- 
cording to each user’s position, comprising: 
a global network of a plurality of reference stations 
for tracking a plurality of satellites (SV’s) in the 
GPS constellation of satellites wherein at least five 
reference stations track each SV; 
a central processing facility (CPF) having means for 
computing satellite orbit and clock errors; and 
communication means connected to the CPF for 
obtaining reference station measurements and for 
distributing DGPS corrections to users; 
wherein the CPF comprises means for building SV 
orbits generating real-time ephemerides in a two 
step least-squares estimator which estimates the 
XYZ position of SV’s and then estimates their 

orbital ephemeris terms. 
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